fastq格式文件处理大全（六）

从计算机的角度来说，生物的序列属于一种字符串，也是一种文本，因此生物信息分析属于文本处理范畴。文本存储为固定格式文件，生物信息的工作就是各种文本文件之间格式的转换，例如通过序列拼接将fastq转换为fasta，通过短序列比对将fastq与fasta合并为bam，通过变异检测将bam中突变位点提取出来转换为vcf。因此，我们可以通过总结每一种生物数据文件格式的处理方法来学习生物信息，这样当拿到固定格式的文件之后，就知道该如何来处理了。

提取序列

建立索引之后就可以从fastq文件中快速提取序列了，不过这个操作一般是对fasta格式文件操作，fastq文件用的并不多，如果需要对fastq文件建立索引，可以使用samtools fqidx命令，不过这里需要注意，fastq格式文件需要使用bgzip压缩，一般的fastq都是采用gzip压缩，需要重新处理文件才行。

#建立索引
samtools fqidx SRR8651554_1.fastq
#提取序列，将ID添加到结尾处，这里提取三条
samtools fqidx SRR8651554_1.fastq  SRR8651554.1 SRR8651554.2 SRR8651554.3

合并两条序列

双末端测序的reads来自一条片段的两侧，如果文库大小比较小，测序读长比较长，例如pairend 300，insertsize 500，就有一些片段中间具有overlap区域，因此可以将两条reads进行拼接，连成更长的片段。练成更长的片段有利于进行序列拼接，也具有更好的唯一性，例如16S序列分析中常用此步骤。有很多工具可以完成这项操作，例如cope，flash，fastq-join等。

cope -a SRR8651554_1.fastq.gz -b SRR8651554_2.fastq.gz -o connect

kmer计数

kmer是一段序列，一般将fastq文件按照固定长度（奇数），例如15，从头到尾按照步移数为1开始截取序列，例如1-15,2-16,3-17……然后对kmer频率进行计数，kmer分析可以用来估计基因组大小等，通过kmer频数分布也可以反应测序错误，或者杂合位点的分布。一般认为kmer频数为1的序列包含测序错误位点。可以使用jellyfish对fastq文件进行kmer计数。

gunzip  SRR8651554_1.fastq.gz
jellyfish count -m 15 -s 100M  -o kmer_counts SRR8651554_1.fastq
jellyfish histo kmer_counts

拼接基因组

质控过滤完的fastq文件可以直接使用拼接软件进行基因组拼接，输入文件为fastq格式，拼接完得到基因组文件，为fasta格式。

python spades.py -o result -1 clean.1.fq.gz -2 clean.2.fq.gz >spades.log 2>spades.err

fastq到bam

很多分析的第一步就是将fastq文件转换为bam，包括变异检测，RNAseq等，如何将fastq文件转换为bam呢，这就需要通过短序列比对。一些测序仪直接输出bam格式文件，例如Ion Torrent，其实那个并不是比对之后得到的bam，其实属于uBam格式。可以将fastq进行短序列比对的工具很多，这里以bwa为例。除了需要fastq格式，还需要fasta格式作为参考序列。

 #建立索引  
bwa index -a is  ref.fna 
#bwa比对
bwa mem -t 4 -R '@RG\tID:A1\tPL:illumina\tSM:MTB' ref.fna  clean.1.fq.gz  clean.2.fq.gz  >A1.sam