诺禾致源有参转录组分析流程资料下载.pdf

1、FASTQ格式文件中每个read由四行描述，如下：HWI-ST1276:71:C1162ACXX:1:1101:1208:2458 1:N:0:CGATGT NAAGAACACGTTCGGTCACCTCAGCACACTTGTGAATGTCATGGGATCCAT+#55?BBBBB?BADEEFFCFFHHFFCFFHHHHHHHFAE0ECFFD/AEHH其中第一行以“”开头，随后为Illumina 测序标识符（Sequence Identifiers）和描述文字（选择性部分）；第二行是碱基序列；第三行以“+”开头，随后为Illumina 测序标识符（选择性部分）；第四行是对应碱基的测序质量，

2、该行中每个字符对应的 ASCII 值减去 33，即为对应第二行碱基的测序质量值。Illumina测序标识符详细信息如下：HWI-ST1276Instrument unique identifier of the sequencer71run number Run number on instrumentC1162ACXXFlowCell ID ID of flowcell1LaneNumber positive integer1101TileNumber positive integer1208X x coordinate of the spot.Integer which can be ne

3、gative2458Y y coordinate of the spot.Integer which can be negative1ReadNumber-1 for single reads;1 or 2 for paired endsNwhether it is filtered-NB：Y if the read is filtered out,not in the delivered fastq file,N otherwise0control number-0 when none of the control bits are on,otherwise it is an even nu

4、mberCGATGTIllumina index sequences5/49北京诺禾致源生物信息科技有限公司2测序数据质量评估测序数据质量评估2.1测序错误率分布检查测序错误率分布检查每个碱基测序错误率是通过测序Phred数值（Phred score,Qphred）通过公式（公式1：Qphred=-10log10（e）转化得到，而Phred 数值是在碱基识别（Base Calling）过程中通过一种概率模型计算得到，这种模型可以准确地预测碱基判别的错误率。Phred分值，不正确的碱基识别率，碱基正确识别率以及Q-score的对应关系如下表所显示：illumina Casava 1.8版本碱基识

5、别与Phred分值之间的简明对应关系Phred分值分值不正确的碱基识别不正确的碱基识别碱基正确识别率碱基正确识别率Q-sorce101/1090%Q10201/10099%Q20301/100099.9%Q30401/1000099.99%Q40测序错误率与碱基质量有关，受测序仪本身、测序试剂、样品等多个因素共同影响。对于RNA-seq技术，测序错误率分布具有两个特点：（1）测序错误率会随着测序序列（Sequenced Reads）长度的增加而升高，这是由于测序过程中化学试剂的消耗导致的，并且为illumina高通量测序平台都具有的特征。（2）前6个碱基的位置也会发生较高的测序错误率，而这个长

6、度也正好等于在RNA-seq建库过程中反转录所需要的随机引物的长度。所以前6个碱基测序错误率较高的原因为随机引物和RNA模版的不完全结合（Jiang et al.）。图1测序错误率分布图横坐标为reads的碱基位置，纵坐标为单碱基错误率6/49北京诺禾致源生物信息科技有限公司2.2A/T/G/C 含量分布检查含量分布检查GC含量分布检查用于检测有无AT、GC 分离现象，而这种现象可能是测序或者建库所带来的，并且会影响后续的定量分析。在illumina测序平台的转录组测序中，反转录成cDNA时所用的6bp 的随机引物会引起前几个位置的核苷酸组成存在一定的偏好性。而这种偏好性与测序的物种和实验室环

7、境无关，但会影响转录组测序的均一性程度（Hansen et al.）。除此之外，理论上普通文库的G和C碱基及A和T碱基含量每个测序循环上应分别相等，且整个测序过程稳定不变，呈水平线，而对于链特异性建库会出现GC分离的现象。对于DGE测序来说，由于随机引物扩增偏差等原因，常常会导致在测序得到的每个read前6-7个碱基有较大的波动，这种波动属于正常情况。图2GC含量分布图横坐标为reads的碱基位置，纵坐标为单碱基所占的比例；不同颜色代表不同的碱基类型7/49北京诺禾致源生物信息科技有限公司2.3测序数据过滤测序数据过滤测序得到的原始测序序列，里面含有带接头的、低质量的reads，为了保证信息分

8、析质量，必须对raw reads进行过滤，得到clean reads，后续分析都基于cleanreads。数据处理的步骤如下：（1）去除带接头（adapter）的reads；（2）去除N（N表示无法确定碱基信息）的比例大于10%的reads；（3）去除低质量reads（质量值 sQ=5 的碱基数占整个 read 长度的 50以上的 reads）。RNA-seq 的接头（Adapter,Oligonucleotide sequences for TruSeqTM RNA and DNA Sample Prep Kits）信息：RNA 5 Adapter（RA5）,part#15013205：5-

9、AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT-3RNA 3 Adapter（RA3）,part#15013207：5-GATCGGAAGAGCACACGTCTGAACTCCAGTCAC（6位index）ATCTCGTATGCCGTCTTCTGCTTG-3图2.3原始数据组成不同颜色的比例分别代表不同成分比例（1）Adapter related：因有接头，过滤掉的 reads数及其占总 raw reads数的比例。（2）Containing N：因 N含量超过 10%，过滤掉的 reads数及其占总raw reads数的比例

10、。（3）Low quality：因低质量，过滤掉的reads数及其占总raw reads数的比例。（4）Clean reads：最终得到的 clean reads 数及其占总 raw reads 数的比例。8/49北京诺禾致源生物信息科技有限公司2.4测序数据质量情况汇总测序数据质量情况汇总样品测序产出数据质量评估情况详见表表1。表1数据产出质量情况一览表Sample nameRaw readsClean readsClean basesError rate（%）Q20（%）Q30（%）GC content（%）sampleA1_134200519334833964.19G0.0395.619

11、1.3248.43sampleA1_234200519334833964.19G0.0494.1389.0748.47sampleA2_134245266331053654.14G0.0395.8891.8748.42sampleA2_234245266331053654.14G0.0493.9888.8448.44sampleB1_132687612313616593.92G0.0395.9091.8849.30sampleB1_232687612313616593.92G0.0494.2389.2449.30sampleB2_130232747292372673.65G0.0395.929

12、1.9248.87sampleB2_230232747292372673.65G0.0494.0588.9448.88sampleC1_128782461283694583.55G0.0395.9591.9848.57sampleC1_228782461283694583.55G0.0494.3989.5048.57sampleC2_128521158280358773.5G0.0395.9492.0047.72sampleC2_228521158280358773.5G0.0494.2889.3847.72数据质量情况详细内容如下：（1）Sample name：样品名，1为左端reads，2

13、为右端reads。样品的 clean reads 总数为 左端+右端。（2）Raw reads：统计原始序列数据，以四行为一个单位，统计每个文件的测序序列的个数。（3）Clean reads：计算方法同 Raw Reads，只是统计的文件为过滤后的测序数据。后续的生物信息分析都是基于Clean reads。（4）Clean bases：Clean reads的个数乘以长度，并转化为以G为单位。（5）Error rate：通过公式1计算得到。（6）Q20、Q30：分别计算 Phred 数值大于20、30的碱基占总体碱基的百分比。（7）GC content：计算碱基G和C的数量总和占总的碱基数量的

14、百分比。9/49北京诺禾致源生物信息科技有限公司2.5 质量评估质量评估Q&A问问：测序错误率会随着测序序列长度的增加而升高，错误率在多少是可以接受的范围？答答：诺禾的测序会进行严格的数据质量把控。一般情况下，单个碱基位置的测序错误率应该低于1%，最高在6%左右可以接受。问问：诺禾质控的标准是什么？是否严格？为保证后续分析的质量，诺禾会严格把控cleandata的筛选标准，具体标准如下：（3）去除低质量reads（质量值sQ 达水平。Reads计数除了与基因的真实表达水平成正比外，还与基因的长度和测序深度成正相关。为了使不同基因、不同实验间估计的基因表达水平具有可比性，人们引入了FPKM的概念

15、，FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairssequenced）是每百万fragments中来自某一基因每千碱基长度的fragments数目，其同时考虑了测序深度和基因长度对fragments计数的影响，是目前最为常用的基因表达水平估算方法（Trapnell,Cole,et al.,2010）。本次采用HTSeq软件对各样品进行基因表达水平分析，使用的模型为union。结果文件分别统计了不同表达水平下基因的数量以及单个基因的表达水平。一般情况下，FPK

16、M数值0.1或者1作为判断基因是否表达的阈值，不同的文献所采用的阈值不同。表7.1不同表达水平区间的基因数量统计表FPKM IntervalsampleA1sampleA2sampleB1sampleB2sampleC1sampleC20110422（43.42%）10393（43.30%）10549（43.95%）10590（44.12%）10423（43.43%）10373（43.22%）132210（9.21%）2203（9.18%）2204（9.18%）2197（9.15%）2267（9.45%）2244（9.35%）3154915（20.48%）4980（20.75%）4926（20

17、.52%）4961（20.67%）4907（20.45%）5084（21.18%）15604586（19.11%）4546（18.94%）4501（18.75%）4466（18.61%）4577（19.07%）4541（18.92%）601867（7.78%）1878（7.83%）1820（7.58%）1786（7.44%）1826（7.61%）1758（7.33%）表7.2基因表达水平统计表Gene_idsampleA1sampleA2sampleB1sampleB2sampleC1sampleC2ENSMUSG00000095309000000ENSMUSG00000029064566.3

18、10996987737575.11185780644426.623104830895427.399760864362.671217091939387.266830611945Novel014397.008204950613217.499303187195.832011182695646.865009730575529.5594181376769911.5285845636247ENSMUSG000000233482.133385703326562.370301894497952.083598558961582.295586491128513.191627149169412.9613042332

19、13426/49北京诺禾致源生物信息科技有限公司7.2 基因表达水平分析基因表达水平分析Q&基因表达水平如何计算？在RNA-seq技术中，FPKM（expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced）是每百万fragments中来自某一基因每千碱基长度的fragments数目,FPKM 时考虑了测序深度和基因长度对reads计数的影响，是目前最为常用的基因表达水平估算方法。认为基因表达的阈值是多少？为什么设置为这个阈值？有参转录组当中，认为FPKM1

20、是基因表达的。这个阈值是主流杂志推荐的，也能够很好的反应基因的表达水平。基因表达水平分析相关名词的解释基因表达水平分析相关名词的解释:FPKM：expected number of Fragments Per Kilobase of transcript sequence per Millions base pairs sequenced，是每百万fragments中来自某一基因每千碱基长度的fragments数目。RPKM：expected number of reads Per Kilobase of transcript sequence per Millions base pairs

21、sequenced，是每百万reads中来自某一基因每千碱基长度的reads数目。27/49北京诺禾致源生物信息科技有限公司8RNA-seq整体质量评估整体质量评估8.1 RNA-Seq相关性检查相关性检查生物学重复是任何生物学实验所必须的，高通量测序技术也不例外（Hansen et al.）。生物学重复主要有两个用途：一个是证明所涉及的生物学实验操作是可以重复的且变异不大，另一个是为了确保后续的差异基因分析得到更可靠的结果。样品间基因表达水平相关性是检验实验可靠性和样本选择是否合理的重要指标。相关系数越接近1，表明样品之间表达模式的相似度越高。Encode计划建议皮尔逊相关系数的平方（R2）

22、大于0.92（理想的取样和实验条件下）。具体的项目操作中，我们要求R2至少要大于0.8，否则需要对样品做出合适的解释，或者重新进行实验。图8.1RNA-Seq相关性检查热图：样品间相关系数热图；散点图（若样品多于4组，则仅展示生物学重复之间的散点图）：样品间的相关系数散点图，R2:pearson相关系数的平方。28/49北京诺禾致源生物信息科技有限公司8.2 RNA-seq整体质量评估整体质量评估Q&样品间的相关性有何意义？如何计算?样品间的相关性反应了样品间的相似程度，即不同处理或组织的样品在表达水平方面的相似度。相关系数越接近1，样品间的相似度越高，样品间的差异基因也越少。生物学重复间的样

23、品的相关系数应大于生物学重复外的样品的相关系数。相关系数的计算方法有三种：A.Pearson correlation;B.Spearman rankcorrelation;C.Kendalls。诺禾使用R语言进行Pearson相关系数的计算。29/49北京诺禾致源生物信息科技有限公司9 差异表达分析差异表达分析9.1基因表达水平对比基因表达水平对比通过所有基因的FPKM分布图以及盒形图对不同实验条件下的基因表达水平进行比较。对于同一实验条件下的重复样品，最终的FPKM为所有重复数据的平均值。图9.1不同实验条件下基因表达水平比对图图一：FPKM盒形图，横坐标为样品名称，纵坐标为log10（FP

24、KM+1），每个区域的盒形图对五个统计量（至上而下分别为最大值，上四分位数，中值，下四分位数和最小值）图二：FPKM分布图，横坐标为log10（FPKM+1）,纵坐标为基因的密度30/49北京诺禾致源生物信息科技有限公司9.2差异表达基因列表差异表达基因列表基因差异表达的输入数据为基因表达水平分析中得到的readcount数据。对于有生物学重复的样品，我们采用DESeq（Anders et al,2010）进行分析：该分析方法基于的模型是负二项分布，第 i 个基因在第 j 个样本中的 read count 值为Kij，则有Kij NB（ij,ij2）对于无生物学重复的样品，先采用TMM对rea

25、d count数据进行标准化处理，之后用DEGseq进行差异分析。差异表达基因列表如下：表9.2差异基因列表Gene IdsampleAsampleBlog2FoldChangepvalp-adjustedENSMUSG0000002450127459.76672652123588.044127161872.9361.5913e-1703.2181e-166ENSMUSG0000003142515394.7853713593106400.873598884-2.7899.2315e-1086.223e-104ENSMUSG00000076439497.8389515519293633.6179

26、2454481-2.86776.3681e-1086.223e-104ENSMUSG0000002286510422.51757993082060.393732559542.33879.3935e-1064.7491e-102差异基因列表主要包括的内容：（1）Gene id:基因编号（2）Sample1：校正后样品1的readcount值（3）Sample2：校正后样品2的readcount值（4）log2FoldChange:log2（Sample1/Sample2）（5）pvalue（pval）:统计学差异显著性检验指标（6）qvalue（p-adjusted）:校正后的pvalue。qv

27、alue越小，表示基因表达差异越显著31/49北京诺禾致源生物信息科技有限公司9.3差异表达基因筛选差异表达基因筛选用火山图可以推断差异基因的整体分布情况，对于无生物学重复的实验，为消除生物学变异，从差异倍数和显著水平两个方面进行评估，对差异基因进行筛选，阈值设定一般为:|log2（FoldChange）|1 且 qvalue 0.005。对于有生物学重复的实验，由于DESeq已经进行了生物学变异的消除，我们对差异基因筛选的标准一般为:padj 0.05。图9.3差异基因火山图有显著性差异表达的基因用红色点（上调）和绿色点（下调）表示，无显著性差异表达的基因用蓝色点表示；横坐标代表基因在不同样

28、本中表达倍数变化；纵坐标代表基因表达量变化差异的统计学显著性32/49北京诺禾致源生物信息科技有限公司9.4差异基因聚类分析差异基因聚类分析聚类分析用于判断差异基因在不同实验条件下的表达模式；通过将表达模式相同或相近的基因聚集成类，从而识别未知基因的功能或已知基因的未知功能；因为这些同类的基因可能具有相似的功能，或是共同参与同一代谢过程或细胞通路。以不同实验条件下的差异基因的FPKM值为表达水平，做层次聚类（hierarchicalclustering）分析，不同颜色的区域代表不同的聚类分组信息，同组内的基因表达模式相近，可能具有相似的功能或参与相同的生物学过程。除了差异基因表达量FPKM层次聚类分析，我们还分别用H-clu