生物数据格式 - fastq

格式

fastq格式是一种包含质量值的序列文件，一般用来存储原始测序文件，文件扩展名一般为fastq或fq，目前主流测序仪器都以fastq格式存储测序数据。fastq的序列格式如下，每条序列的信息包括四行。 fastq文件格式.png

第一行：以’@‘开头，是这条read的名字，它是每一条read的唯一标识符，在同一份fastq文件中不会重复出现，但在双端测序中，一条read的名字在R1和R2中可能都存在；
第二行：有A、G、C、T和N五种字母构成，是测序read的序列，其中N代表的是测序时无法被识别的碱基；
第三行：以‘+’开头，备注行，在旧版fastq文件中会直接重复第一行的信息，现在为节省存储空间一般不加；
第四行：测序序列的质量信息，俗称“Q值”，指每一个碱基的可靠程度，用ASCII码表示。Q值的计算公式：Q=-10lgP，P:碱基测错概率。Q值加上33后转换成ASCII码（需要一个Q值符号对应一个碱基；33以前的ASCII码无法用键盘字符表示；illumina测序1.8版本以上加33，以下加64。）fastq格式中质量值具有重要的作用，在很多的分析例如数据质控，数据过滤，序列拼接，短序列比对，变异检测中都会被用到。

获取

除了直接对核酸测序获得fastq文件以外，还可以在NCBI SRA数据库下载开放的测序文件，通过sratools工具中的prefetch可以下载fastq文件

# SRA下载数据
## prefetch/wget下载.sra文件
prefetch SRR号
wget 数据超链接地址
## fastq-dump 将sra文件转换成fastq格式
fastq-dump --split-3 --gzip -A 样本名 XX.sra #--split-3参数可以讲sra文件拆分成双端的fastq，--gzip参数设置压缩生成的文件，-A参数设置样本名

# 直接测序数据
## 采用gzip压缩或解压缩
gzip XX.fastq #压缩
gzip -d XX.fastq.gz #解压缩
md5sum -c SRR8651554_1.fastq.md5 #通过md5值校验文件完整性，防止文件传输不完整

处理

处理文件

## 合并
cat 1.fastq 2.fastq > new.fastq #对象未压缩，末尾追加
zcat 1.fastq.gz 2.fastq.gz > new.fastq.gz #压缩对象，末尾追加
seqtk mergepe 1.fastq.gz 2.fastq.gz > new.fastq.gz #压缩对象，逐条合并
## 拆分 使用seqkit split按ID拆分或者指定切割份数或大小
seqkit split [flags] #Flags：-i,按ID切割；-p,切割成N份；-s,按N条每份切割
## 排序 使用seqkit sort通过id/名称/序列/长度来排序
seqkit sort [flags] #Flags: -l,通过序列长度；-n,通过全名而不是id；-s,通过序列；-i,忽视大小写；-r,反向排序；-2,双流程模式读取文件来降低内存使用。
## 抽样 
seqtk seq -f 0.1 -s 11 XX.fastq.gz #seqtk 按照百分比
seqkit sample -p 0.1 XX.fastq.gz #seqkit 按照百分比
## 提取 建立索引后根据ID快速提取序列，适用fasta文件
samtools fqidx XX.fastq #建立索引
samtools fqidx XX.fastq  ID1 ID2 ID3 > new.fastq #提取序列，
## 转fasta
seqtk seq -A XX.fastq.gz #seqtk工具
seqkit fq2fa XX.fastq.gz #seqkit工具
## 列表转换
seqkit fx2tab XX.fastq.gz #将fastq格式转换为列表格式，将四行数据转换为一行三列，方便根据ID进行处理。
seqkit tab2fx XX.table #将列表转为换fastq格式
## 质量转换 之前测序的文件可能是phred+64的模式，可以转换成现在phred+33的格式。
seqkit convert --to Illumina-1.5+ XX.fastq.gz |head -4 #将illumina 1.8转换为1.5
seqkit convert  XX_illmina1.5.gz #将illumina 1.5转换为1.8，什么都不加就是转换为1.8

处理序列

## 统计
seqkit stats XX1.fastq.gz XX2.fastq.gz #统计序列条数，碱基总数，reads读长分布等
seqtk comp XX1.fastq.gz XX2.fastq.gz #统计fastq文件每条序列ATCG四种碱基组成以及质量值分布
seqtk fqchk XX1.fastq.gz #获取每个碱基点的分布和质量值，结果可以用来绘制碱基质量以及含量分布图
## 质控QC
fastqc -f fasqc *.fq.gz #-f参数可以设置对fastq、bam、sam文件进行质控
## 过滤短序列 Ion Torrent，pacbio，nanopore测序的fastq文件序列长度并不相同，通常需要过滤较短的序列。
seqkit seq -m 150 XX.fastq.gz | gzip -  >XX_filter_150.fq.gz #过滤小于150bp序列，并压缩输出 
seqtk seq -L 150 XX.fastq.gz #过滤小于150bp序列
seqkit seq -M 150 XX.fastq.gz #保留小于150bp序列
## 去接头
cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC test.fq > output.fq #单端测序，3‘端精确匹配上a后序列的read被去除
cutadapt -a GATCGGAAGAGCACACGTCTGAACTCCAGTCAC -O 5 test.fq -o test_output.fq --discard-trimmed #单端测序，3‘端匹配上a后序列5个碱基的read被去除
cutadapt -a GGAAGGTGGGGATGACGT -A AACMGGATTAGATACCCKG -o output1.fastq -p output2.fastq input_1.fastq.gz input_2.fastq.gz #双端测序，参数：-a,read1的adapter; -A,read2的adapter; -o,read1的输出; -p,read2的输出
## 去低质量
seqtk seq -q 20 XX.fastq.gz |less -S #过滤含质量q值小于20的reads
## 截取头尾
seqtk trimfq -b 15 -e 15 -Q XX.fastq.gz  #参数：-b,开头; -e,尾部
## 过滤 
fastp -i XX_1.fastq.gz -I XX_2.fastq.gz -o clean.1.fq.gz -O clean.2.fq.gz -z 4 -q 20 -u 40 -n 10 -f 15 -t 15 -F 15 -T 15 -h fastp.html 
 #   -z 5 设置gzip压缩级别，默认4,1~9 1快-9慢
 #   -q, --qualified_quality_phred 设置碱基质量值不小于多少时，该碱基为合格碱基，默认碱基质量值是15，即默认碱基质量>=15是合格碱基，<15为不合格碱基；
 #   -u, --unqualified_percent_limit 设置允许不合格碱基的占比为多少时，去掉这条read，默认是40，即默认不合格碱基占比>40%时，去掉该read；
 #   -n 过滤 N碱基过多的reads，如果N碱基含量大于n，这条read/pair将被舍弃，默认5；
 #   -f, --trim_front1 对R1起始几bp进行去除，例如：-f 15或--trim_front1=15表示去除R1起始位置15bp碱基；
 #   -t, --trim_tail1 对R1末尾几bp进行去除，例如：-t 15或--trim_tail1=15表示去除R1末尾位置15bp碱基；
 #   -F, --trim_front2 与R1相似；不设置默认则与R1指定的参数相同；
 #   -T, --trim_tail2 与R1相似；不设置默认则与R1指定的参数相同；
 #   -h, --html 设置输出html格式的质控结果文件名，不设置则默认html文件名为fastp.html
## 拼接
### 某些大小比较小，测序读长比较长的文库有一些片段的两条reads中间具有overlap区域，因此可以将两条reads进行拼接，连成更长的片段，这样有利于进行序列拼接，也具有更好的唯一性。
cope -a XX_1.fastq.gz -b XX_2.fastq.gz -o connect
## kmer计数
### kmer是一段序列，一般将fastq文件按照固定长度（奇数），例如15，从头到尾按照步移数为1开始截取序列，例如1-15,2-16,3-17……然后对kmer频率进行计数，kmer分析可以用来估计基因组大小等，通过kmer频数分布也可以反应测序错误，或者杂合位点的分布。一般认为kmer频数为1的序列包含测序错误位点。可以使用jellyfish对fastq文件进行kmer计数。
jellyfish count -m 15 -s 100M  -o kmer_counts XX.fastq
jellyfish histo kmer_counts
## 拼接基因组
### 质控过滤完的fastq文件可以直接使用软件进行基因组拼接，输入文件为fastq格式，输出为fasta格式。
python spades.py -o result -1 XX_clean.1.fq.gz -2 XX_clean.2.fq.gz >spades.log 2>spades.err
## 转bam
### 包括变异检测，RNAseq等很多分析的第一步就是将fastq文件转换为bam，这就需要通过短序列比对，有多种软件可以完成这个过程。(一些测序仪直接输出bam格式文件，例如Ion Torrent，其实那个并不是比对之后得到的bam，其实属于uBam格式。)
bwa index -a is  ref.fa  #建立索引  
bwa mem -t 4 -R '@RG\tID:A1\tPL:illumina\tSM:MTB' ref.fa  clean.1.fq.gz  clean.2.fq.gz  >XX.sam #bwa比对 
    #   -t int  线程数，默认是1。
    #   -R str  完整的read group的头部，可以用 '\t' 作为分隔符，在输出的SAM文件中被解释为制表符TAB. read group 的ID，会被添加到输出文件的每一个read的头部。

———以上纯属个人理解与记录