分类： 生物信息

fastaq文件中reads数的统计

在汇总测序的基本信息时经常需要从fastaq文件中计算某个样本的原始的read数量。

fastaq文件的一般格式为：

@FCD056DACXX:3:1101:2163:1959#TCGCCGTG/1
TCCGATAACGCTCAACCAGAGGGCTGCCAGCTCCGATCGGCAGTTGCAACCCATTGGCCGTCTGAGCCAGCAACCCCGGA
+
gggiiiiiiiiiiiiiiiiiiiiiiiiiigggggeeecccccc^bcbcccccccbccccc]aaccbbccc^R^^acccc_

每四行一条记录， 每条记录以@符号开始，后面跟着的是该reads的标识符，第二行为对应的核酸序列，第三行为‘+’，第四行为序列的质量信息，以字符对应的ASCII码来表示对应碱基的测序质量。因此，想要知道有多少个reads只要统计有多少个只包含+的行就知道有多少个reads了。如下：

cat your.fastq.gz | grep -c '^+$'

有些人给出的方法是：cat your.fastq.gz | grep -c ‘+’  或 cat your.fastq.gz | grep -c ‘^+’ 这两种方法都是不严谨的，因为第四行的质量信息中也是包含了‘+’的，如果只是统计包含‘+’或以‘+’开头的行会多统计出一部分数据来，尤其是第一种情况。因此，正确的方式应该是计数只包含‘+’的行。此外要注意的是，如果是pair-end的数据，一个样本是包含了两个文件的，要放在一起计数。这和后面统计质控后的Clean reads以及计算后面比对reads的比例有关系。

测序基本信息：比对深度和覆盖度

测序的比对率、测序深度和覆盖度是了解测序情况的最基本信息。很多工具都能完成这些数据的分析。现在对网上常用的方式收集汇总了一下，按照软件一个一个滴说。

1. Samtools

#命令解析 
$ samtools flagstat sample1.bam 
3826122 + 0 in total (QC-passed reads + QC-failed reads) #总共的reads数 
0 + 0 secondary 
1658 + 0 supplementary 
343028 + 0 duplicates 
3824649 + 0 mapped (99.96% : N/A) #总体reads的匹配率 
3824464 + 0 paired in sequencing #总共的reads数 
1912442 + 0 read1 #reads1中的reads数 
1912022 + 0 read2 #reads2中的reads数 
3808606 + 0 properly paired (99.59% : N/A) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值 
3821518 + 0 with itself and mate mapped #paired reads中两条都比对到参考序列上的reads数 
1473 + 0 singletons (0.04% : N/A) #单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。 
5882 + 0 with mate mapped to a different chr#paired reads中两条分别比对到两条不同的参考序列的reads数 
4273 + 0 with mate mapped to a different chr (mapQ>=5) #paired reads中两条分别比对到两条不同的参考序列的reads数

samtools常用命令详解

转载自：https://www.cnblogs.com/emanlee/p/4316581.html

samtools的说明文档：http://samtools.sourceforge.net/samtools.shtml
samtools是一个用于操作sam和bam文件的工具合集。包含有许多命令。以下是常用命令的介绍

1. view

view命令的主要功能是：将sam文件转换成bam文件；然后对bam文件进行各种操作，比如数据的排序(不属于本命令的功能)和提取(这些操作 是对bam文件进行的，因而当输入为sam文件的时候，不能进行该操作)；最后将排序或提取得到的数据输出为bam或sam（默认的）格式。

bam文件优点：bam文件为二进制文件，占用的磁盘空间比sam文本文件小；利用bam二进制文件的运算速度快。

view命令中，对sam文件头部的输入(-t或-T）和输出(-h)是单独的一些参数来控制的。

Usage: samtools view [options] | [region1 [...]]
默认情况下不加 region，则是输出所有的 region.

Options: -b       output BAM
                  默认下输出是 SAM 格式文件，该参数设置输出 BAM 格式
         -h       print header for the SAM output
                  默认下输出的 sam 格式文件不带 header，该参数设定输出sam文件时带 header 信息
         -H       print header only (no alignments)
         -S       input is SAM
                  默认下输入是 BAM 文件，若是输入是 SAM 文件，则最好加该参数，否则有时候会报错。
         -u       uncompressed BAM output (force -b)
                  该参数的使用需要有-b参数，能节约时间，但是需要更多磁盘空间。
         -c       Instead of printing the alignments, only count them and print the 
                  total number. All filter options, such as ‘-f’, ‘-F’ and ‘-q’ , 
                  are taken into account.
         -1       fast compression (force -b)
         -x       output FLAG in HEX (samtools-C specific)
         -X       output FLAG in string (samtools-C specific)
         -c       print only the count of matching records
         -L FILE  output alignments overlapping the input BED FILE [null]
         -t FILE  list of reference names and lengths (force -S) [null]
                  使用一个list文件来作为header的输入
         -T FILE  reference sequence file (force -S) [null]
                  使用序列fasta文件作为header的输入
         -o FILE  output file name [stdout]
         -R FILE  list of read groups to be outputted [null]
         -f INT   required flag, 0 for unset [0]
         -F INT   filtering flag, 0 for unset [0] 
                  Skip alignments with bits present in INT [0]
                  数字4代表该序列没有比对到参考序列上
                  数字8代表该序列的mate序列没有比对到参考序列上
         -q INT   minimum mapping quality [0]
         -l STR   only output reads in library STR [null]
         -r STR   only output reads in read group STR [null]
         -s FLOAT fraction of templates to subsample; integer part as seed [-1]
         -?       longer help

继续阅读“samtools常用命令详解”