一、测序数据

1.Raw data（原始数据）

公司一次测序产生的全部原始数据。理论上，它们应该是没有经过任何过滤的，无论好坏。

image.png

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和后续的基础分析结果。

二、测序数据信息的统计

image.png

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)