linux质控命令,Biostar_handbook||Charpter_6789_数据的格式_获取_质控_seqkit

Charpter6：数据的格式Data Formats

常用数据库

NCBI

EBI

Uniprot

Phytozome

Ensemble plants

TAIR

PlantGDB

PlantTFDB 植物转录因子数据库

这些是我目前常用的一些数据库，可以看出来我是做植物的。其实植物领域还有许多重要的数据库，更不必说动物的了，各个物种都能拿出来做个数据库。数据这么多，如何更加高效的从数据库中实现你想要的目的是门学问。而自己目前对这些数据库知道的也只是初步的一些东西。其实除了NCBI外很多重要的数据库都有一些很实用的功能，但是却缺少教程指导推广，以后有机会可以再对一些数据库的功能进行探索探索。

数据格式

1. GeneBank format

最常见的生物数据存储方式FASTA,GeneBank, FASTQ,前两者被称为curated sequence information.

Genebank/gb格式文件为NCBI存储基因信息的格式。其中NCBI里的RefSeq数据库为标准非冗余的基因数据库.(It provides the Sequences Record and Baseline for biological studies and it's a non-redundant set of reference standards derived)

用readseq对gb数据转换：

efetch -db=nuccore -format=gb -id=AF086833 > AF086833.gb

###转gb到cds序列fasta序列

readseq -p -format=FASTA AF086833.gb ## complete genome seq

readseq -p -format=FASTA -feat=CDS AF086833.gb

readseq -p -format=GFF AF086833.gb

readseq -p -format=GFF -field=CDS AF086833.gb

2. FASTA

就是第一行以'>'开始的序列的各种信息，第二行具体atgc

3. FASTQ

二代测序数据存储格式，第一行'@'，第二行''序列，第三行'+或者序列信息''，第四行测序的数据质量。

Charpter 7: 获取序列How to get DATA

1. 通过LINUX系统直接从NCBI上下载数据

### 以genebank格式序列下载

efetch -db=nuccore -format=gb -id=AF086833 > AF086833.gb

### 以fatst格式下载序列

efetch -db=nuccore -format=fasta -id=AF086833 > AF086833.fasta

### 取序列的某一段

efetch -db=nuccore -format=fasta -id=AF086833 -seq_start=1 -seq_stop=100

## 互补序列

-strand 1

-strand 2

通过对文章里的project accession 号来下载数据，以PRJNA257197为例子**

esearch -db nucleotide -query PRJNA257197

esearch -db nucleotide -query PRJNA257197 | efetch -format=fasta > genomes.fa

-db=protein

2. 下载SRA数据库数据

SRA数据库(short read Archive):数据存储的格式

Bioproject：PRJN开头，研究项目的介绍

BioSample：SRS/SRMN开头，生物样品来源的介绍

SRA Experiment: 单独样本的测序文库

SRR Run：直接相关数据存储格式

==根据eseach 和 efetch直接下载数据== (都不用各种网站点了，实用下载数据新技能！！！)

esearch -db sra -query SRP020911 | efetch -format runinfo >SRP020911.csv

## 查看csv第一行标题信息

head -1 SRP020911.csv |tr ',' '\n'|nl ### 第一列为 SRR信息号，第十列为SRR的下载地址

### 获取SRR的信息号,SRR下载地址

cat SRP020911.csv| tail -n +2 | cut -d, -f 1

tail -n +2 | cut -d, -f 10

### 结合xargs直接完成数据的下载：

cat SRP020911.csv| tail -n +2 | cut -d, -f 1 | xargs -i prefetch {}

cat SRP020911.csv| tail -n +2 | cut -d, -f 10 | xargs -i echo wget {} \&

3. 对Seqkit工具的应用

对fq文件基本情况统计seqkit stat *.gz

统计每个序列的GC情况 seqkit fx2tab --name --only-id --gc *.gz

==根据id号提取序列==

### 随机创造个id列表

seqkit sample -p 0.001 duplicated-reads.fq.gz |seqkit seq --name --only-id > id.txt

###根据id列表提取序列

seqkit grep --pattern-file id.txt duplicated-reads.fq.gz

或者-f参数

### 提取目标序列

seqkit grep -p "AT1G30490.1" ath.pep.fa

去除掉含有简并碱基的序列如N/K等

### 生成含有简并碱基序列的id

seqkit fx2tab -n -i -a viral.1.1.genomic.fna.gz | csvtk -H -t grep -f 4 -r -i -p "[^ACGT]" | csvtk -H -t cut -f 1 > id.txt

### 去除id.txt，即反向匹配。

seqkit grep --pattern-file id.txt --invert-match viral.1.1.genomic.fna.gz > clean.fa

### 定位这些含有简并碱基的序列。

seqkit grep -f id2.txt viral.1.1.genomic.fna.gz |seqkit locate -i -P --pattern K+ --pattern N+

定位motif

## 以文件中的motif为标准定位

seqkit locate --degenerate --ignore-case --pattern-file enzymes.fa viral.1.1.genomic.fna.gz

### 根据自己的指定的motif序列定位，可包含简并碱基

seqkit locate -i -d -p CHWWWWWWDG seq.fa

根据序列的长短从小到大排序

seqkit sort --by-length virus.fna.gz > virus.genomics.sorted.fa

文件过大用 --two-pass

Charpter 8：二代三代测序数据基本信息Sequencing instruments

DNA sequencing at 40: past, present and future

1. 测序仪器

一代测序：Sanger法，利用加不同荧光的双脱氧核苷酸ddNTP会中止PCR。

速度快，一次只能一条，1000bp~1500bp左右

二代测序：以Illumina为代表边合成边测序。对随机打断成150~300bp的短片段 +上固定碱基的接头(adapter)+标签(tag/index)。

由于建库中利用了PCR富集序列，因此有一些量少的序列无法被大量扩增，造成一些信息的丢失，且PCR中有概率会引入错配的碱基

三代测序：

单分子荧光测序(SMRT):

纳米孔测序

无需扩增，基于纳米科技，对单分子链DNA/RNA直接用合成/降解/通过纳米孔等方式直接测序

2. 错误率

ILLUMINA 0.1%的错误率，理论上一条序列的测序准确度是递减的，也就是前面准确度高，后面低。且在一些GC率高的区域，测序质量会显著降低。

PacBio 10%的错误率

MinION ：最高达到20%的错误率，且错误多为固定错误，无法通过测序深度来矫正。

可以看到三代测序的数据错误率都很高，所以常需要加测二代的数据来对三代数据矫正。

Charpter 9：测序质控

分为两步：第一步FastQC/MultiQC展示测序的质量；第二步Trimmomatic/fastp去除低质量序列

1. FastQC/MultiQC 显示测序质量

### 下载数据

wget http://data.biostarhandbook.com/data/sequencing-platform-data.tar.gz

tar xzvf sequencing-platform-data.tar.gz

### 直接fastqc显示测序质量

ls *.fq.gz |parallel fastqc -o ./ --nogroup {} &

### MultiQC批量显示测序的质量

##在已经用fastqc得到测序的html文件和zip文件夹

multiqc *fastqc.zip --pdf

在运行过程中曾出现too many levels 错误，需要调用fastqc软件的绝对地址来运行程序

2. Trimmomatic/fastp质量控制

用Trimmomatic 从5‘端开始切除低质量碱基(LEADING：5)TRAILING:20表示切除碱基质量小于20的碱基，MINLEN:50表示过滤切除后长度小于50个碱基的reads

滑动窗口：SLIDINGWINDOW:3：20表示设置窗口大小为3个碱基，从5'端开始切除reads，read的平均质量低于20，开始切除后边的碱基序列。过滤长度小于50碱基的reads，使用更大的窗口可以减轻切除的力度

trimmomatic PE SRR11111_1.fq SRR11111_2.fq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq ILLUMINACLIP:/home/wangtianpeng/anaconda3/share/trimmomatic-0.36-5/adapters/TruSeq3-PE-2.fa:2:30:10 LEADING:5 TRAILING:5 SLIDINGWINDOW:4:5 MINLEN:20

trimmomatic PE SRR1553607_1.fastq SRR1553607_2.fastq trimmed_1.fq unpaired_1.fq trimmed_2.fq unpaired_2.fq SLIDINGWINDOW:4:30

去除特定的adapter

cat ~/miniconda3/envs/bioinfo/opt/fastqc-0.11.5/Configuration/adapter_list.txt

echo ">nextera" > nextera.fa

echo "CTGTCTCTTATACACATCTCCGAGCCCACGAGAC" >> nextera.fa

## 最后加上ILLUMINACLIP:参数

ILLUMINACLIP:adapter.fa:2:30:5

一键化软件fastp

###一键化输出

fastp -i in.R1.fq -o out.R1.fq -I in.R2.fq -O out.R2.fq

### 默认自动化去接头，如果不需要参数-A

### PE 数据的碱基校正 PE 数据的每一对 read 进行分析，查找它们的 overlap 区间，然后对于 overlap 区间中不一致的碱基，如果发现其中一个质量非常高，而另一个非常低，则可以将非常低质量的碱基改为相应的非常高质量值的碱基值

-c 参数使用

fastp -q 30 -5 -l 100 -i il_1.fq.gz -I il_2.fq.gz -o i1_clean_1.fq -O i1_clean_2.fq

这里标准为：平均质量高于Q30，对5‘端进行低质量碱基删除，保留大于100bp的短读。

fastp -3 -W 4 -M 20 -a AGATCGGAAGAG -i $i.fastq.gz -o ../clean_data/fastp_out_2.0/$i.qc.fq.gz

对3'端低质量碱基删除，-W -M 为sliding window参数 -a 指定接头序列