一、测序数据
1.Raw data（原始数据）
- 公司一次测序产生的全部原始数据。理论上，它们应该是没有经过任何过滤的，无论好坏。
2. PF data（PF数据）
在测序过程中，Illumina内置软件根据每个测序片段（read，通常每个片段长100个碱基）前25个碱基的质量决定该read是保留还是抛弃。如果没有达到质控标准，则该read的全部碱基都被抛弃；达到标准、保留下来的数据叫做PF data。 PF代表pass filtering。
3. Q30 data
Illumina内置软件根据统一设定的标准来评判碱基识别结果的可靠性，为每个碱基给予一个质量评分（QV）。PF data里质量评分>=30分的数据称为Q30 data。 Q30的意思是该碱基的可靠性为99.9%。Q30数据通常占PF数据的80%左右。视样本质量、操作水平、试剂质量、仪器状态的不同，这一比例有很大波动。
4. Clean data
某些实验室根据其自身的判断标准，在PF data的基础上，进一步删除质量不好的reads后得到的数据。常见的删除动作有：去接头、去N含量高的reads、去质量评分低的reads、去掉每个read的最后几个碱基，等等。
Clean data是国内叫法；PF data是来自Illumina的概念，是广为接受的国际通行标准。
PF算法实质上是选取每个测序片段（read）前25个碱基的质量来代表整条片段的质量，从而决定该片段的去留。Illumina之所以这样做，而不是逐个检查整条片段所有碱基的质量，一方面是为了节省电脑资源，不致于花费太多时间进行运算，拖累测序进程，另一方面也是在大量测序数据的统计结果基础上选择的平衡点，只要前25个碱基是正常的，后75个碱基出问题的概率比较小。
一次测序实验完成，测序仪上展示的数据量和%Q30都是以PF数据为基础的。只要对数据质量有足够信心，就不会对PF数据再进行加工，可以直接把PF数据交给客户，进行下游的生物信息学分析。一般公司会提供clean data和后续的基础分析结果。
二、测序数据信息的统计
1. Clean reads (millions)
- 计算Clean reads数
vim readfq
1 #reads number
2 #!/bin/bash
3
4 ls /public/home/thu/3.RNA_seq/3.sra_to_fastq/*.fastq.gz >fq.list
5
6 for i in $(cat fq.list)
7 do
8 i=`basename $i`
9 printf $i "\t" >>fq.reads_num
10 readfq $i >> fq.reads_num #readfq函数
11 done
- 结果-得到reads和碱基数
ZT6_RNASeq_rep1.sra_1.fastqNum reads:36189024 Num Bases: 5168633934
ZT6_RNASeq_rep1.sra_2.fastqNum reads:36189024 Num Bases: 5169091324
ZT6_RNASeq_rep2.sra_1.fastqNum reads:25599326 Num Bases: 3562568383
ZT6_RNASeq_rep2.sra_2.fastqNum reads:25599326 Num Bases: 3563756693
control_RNASeq_rep1_1.fastq.gzNum reads:24119971 Num Bases: 3617995650
control_RNASeq_rep1_2.fastq.gzNum reads:24119971 Num Bases: 3617995650
control_RNASeq_rep2_1.fastq.gzNum reads:24041041 Num Bases: 3606156150
control_RNASeq_rep2_2.fastq.gzNum reads:24041041 Num Bases: 3606156150
2. Total mapping rate
- samtools flagstat 函数输出mapping rate
ls *sorted.bam | while read id ;do (samtools flagstat $id > $(basename $id ".sorted.bam").stat);done
- 结果分析
1 40621782 + 0 in total (QC-passed reads + QC-failed reads) #通过质控的总reads数
2 0 + 0 secondary
3 0 + 0 supplementary
4 0 + 0 duplicates
5 40621782 + 0 mapped (100.00% : N/A) #对比到基因组上的 reads数
6 40621782 + 0 paired in sequencing # paired reads 数目
7 20445256 + 0 read1 # read1 的数量
8 20176526 + 0 read2 # read2 的数量
9 37510016 + 0 properly paired (92.34% : N/A) #完美匹配的reads数：比对到同一条参考序列，并且两条reads之间的距离符合设置的阈值
10 37883263 + 0 with itself and mate mapped #paired reads中两条都比对到参考序列上的reads数
11 2738519 + 0 singletons (6.74% : N/A) #单独一条匹配到参考序列上的reads数，和上一个相加，则是总的匹配上的reads数。
12 10639 + 0 with mate mapped to a different chr #paired reads中两条分别比对到两条不同的染色体的reads数
13 10639 + 0 with mate mapped to a different chr (mapQ>=5) #paired reads中两条分别比对到不同染色体的且比对质量值大于5的数量
- 输出结果
for id in *.stat ;do echo -e $id >>test `sed -n "5p" $id` >> test ;done #整理结果输出第五行
cat test
#得到文件
control_RNASeq_rep1.stat 28100200 + 0 mapped (100.00% : N/A)
control_RNASeq_rep2.stat 28246701 + 0 mapped (100.00% : N/A)
ZT6_RNASeq_rep1.stat 40621782 + 0 mapped (100.00% : N/A)
ZT6_RNASeq_rep2.stat 32144781 + 0 mapped (100.00% : N/A)
参考
- raw data/PF data/Q30 data/clean data的不同 - 记号晴 - 博客园 https://www.cnblogs.com/huangyinger/articles/10232967.html
- samtools常用命令详解 - emanlee - 博客园 https://www.cnblogs.com/emanlee/p/4316581.html
- 测序数据基本信息统计 | reads,coverage,depth - 简书 https://www.jianshu.com/p/82ed6e27f571
- Samtools 手册： http://www.htslib.org/
- 如何估算测序数据量？ - 知乎 https://zhuanlan.zhihu.com/p/40040208
