中医药mrna mirna大鼠23samples测序项目结题报告

1目1目 —项目说 二实验流 三实验结 四分析流 五分析结 Reads质量评 ReadsGC含量评 Reads核酸组成偏评 Reads重复评 差异表达GO功能分 差异表达的KEGGpathway功能分 SNP和INDEL分 SNP和INDEL的检 SNP和INDEL的注 2附录1结果文件存放信 附录2分析软件说 3—项目说明3 32 S1*vsS2*；S1*vsS3*；S1*vsS4*；S2*vsS3*；S2*vsS4*;S3*vsS4*； Rat（RattusnorvegicusRnor Reads保留比例为99%以上；平均长度150bp。 reads比对到参考组上差异筛选针对所有样本采用STAR软件统计的原始序列计数，针对有重复样本的实验设计利用DESeq2软件对不同样本组之间筛选差异表达的已知，针对无重复样本的实验设计利用DESeq软件对不同样本组之间筛选差异表达的已知。GO功能富集分 SNP和INDEL析

将全部/转录本作为背景列表，差异/转录本列表作为从背景列表中筛选出来的候选列表，利用Fisher精确检验计算代表GO功能集在差异/转录本列表中是否显著富积的PPBenjamini&Hochberg将全部/转录本作为背景列表，差异/转录本列表作为从背景列表中筛选出来的候选列表，利用Fisher精确检验计算代表KEGG功能集在差异/转录本列表中是否显著富积的PPBenjamini&Hochberg采用BWA软件对预处理序列与物种的参考组序列进行序列比对，针对比对结果文件采用picard-tools于PCRGTAK3分别进行SNP和INDEL 差异表达

差异基因采用皮尔逊相关系数(Pearsoncorrelationcoefficient)且满足R2>=0.9来预测共表达对从而构建共表达网络。4二实验流程45三实验结果56四分析流程67五分析结果7本实验采用illuminaHiseq平台的双端模式对多个样本进行高通量，根据illumina数据的低质量分数集中于末端的分布特点,利用FastQC软件对预处理数据进行质量控制分析，并用fastqc-mcf进行去除adapter，低质量的read等。1.1RawRawreadsRawAverageCleanCleanreadsCleanAverageCleanClean Raw Rawread Raw Clean 8CleanreadpairsCleanbases

Cleanreads Cleanbases 991.1.2图 readGC图 Reads核酸组成偏评估是基于Reads比对组结果进行分析的，用于检测核酸组可从NVC（Nucleotideversuscycle）分析中检测到。该图用于评估所有Reads在每处位置RNA，Reads的每个位置上的期望出现的情况是A%=C%=G%=T%=25%。Reads重复评估是基于Reads比对组结果进行分析的。Reads重复评估采取两种策ReadReads；第二是基于比对的方式，比对到相同组位置的Reads被认定为重复Reads。对于剪切体Reads，起始位置比对一致且剪切方式一致的Reads也同样被认定为重复Reads。结果输出文件路径 针对每个样本，利用HISAT软件将预处理序列与物种的参考组序列进行序列比对，举例比如人的参考组版本是rn6，已知转录本在组位置信息文件版本是表 Mapped Reads Reads Non-Splice ReadsProper表表 Mapped Uniquely Readsmapto Readsmapto Splice ReadsProper 从Reads比对到参考 组结果统计每条的序列数，统计每一条上序列 2.2.1组位置信息和已知转录本在组位置信息可以把序列分为以下2种方式：read1mappedto'+'strandindicatesparentalgeneon'+'strandread1mappedto'-'strandindicatesparentalgeneon'-'strandread2mappedto'+'strandindicatesparentalgeneon'-'strandread2mappedto'-'strandindicatesparentalgeneon'+'strand2）read1mappedto'+'strandindicatesparentalgeneon'-'strandread1mappedto'-'strandindicatesparentalgeneon'+'strandread2mappedto'+'strandindicatesparentalgeneon'+'strandread2mappedto'-'strandindicatesparentalgeneon'-'strandFractionofreadsfailedtodetermine:Fractionofreadsexinedby"1++,1--,2+-,2-+":0.5016Fractionofreadsexinedby"1+-,1-+,2++,2--":0.4949FractionofreadsfailedtodetermineAB的比例分别接近0.5，那么表明实验是RNAseq非链特异性；假如AB的比例其中1个接近于1，那么表明实验是RNAseq链特异性。结果表明该实验是RNAseq非链特异性。计算比对参考组的Reads在不同成分内的覆盖情况(例如CDSexon、5'UTRexon、3'UTRexon、Intron、Intergenicregions)。Reads比对分布评估统计每个区域比对到readsreads子和内含子(两个不同的转录本元件)，则按一定的优先顺序只记录一次该Reads在优先CDSexonsUTRexonsIntronsIntergenicregions排序。例如，如果一个Reads比对到一个属于CDSexon和intron区域，该Reads将被标注为属于CDSexons。表 TotalTagsTotalAssignedTags":序列切割（spliced）一次计数成2个（tags）切割2次计数3个 ，所以"Totaltags">="TotalReads"；"TotalAssignedTags":能正确分配到表格中10组分类的 分配到"TSS_up_1kb"也属于"TSS_up_5kb"和"TSS_up_10kb"，"TotalAssignedTags"=CDS_Exons+5'UTR_Exons+3'UTR_Exons+Introns+TSS_up_10kb+和3UTR；3）TSS10kTES10K 评估方转录本覆盖匀一性分布评估用于检测转录本内reads是否均一并且是否存5/偏差。评估方法是：分析过程中将所有已知转录本归一化为长度100nt区域并计算覆盖在每个子区域位置上的reads数目。最终提供描述体5/方向的reads图 5101595%用于分别计图 可变性剪切位点注释评估分析过程是依据各样本比对参考组与已知模型注释信息，比较已知plceuncton（剪切位点）获得当前转录组内的新pceuncton（剪切位点数目与比例。每个探测到的pceuncton可被划分为3个独立类型：1)Annotat全部属于已知模型注释内的剪切位点、即包括剪切位点的5端剪切位点和3端剪切位点。2)copltenov:全部属于新的剪切位点。剪切位点两端均不属于已知模型中被注释的部分。3)partanove:某部分（5SSor3SS）属于已知模型注释内的剪切位点，另外部分（3SSor5SS）属于新的剪切位点的情形。38353721112531chrom

startpositionofjunction(coordinateis0endpositionofjunction(coordinateis1numberofspliceeventssupportingthis

5101595结果输出文件路径 /针对每个样本，从Reads比对到参考组的结果利用cufflinks软件基于已知转录本在组上位置信息文件作为指导组装样本中的转录本，分别将所有样本的Cufflinks组装好的转录本进行再次组装。组装转录本在组位置信息GTF格式文件见表5.1.1。Cufflinks组装原理如下：首先从比对结果识别不兼容的双端片段(pairsof‘lefragments)，这些片段之间可能是由于不同形式的mRNA剪切异构体(splicedmRNAisoforms)形成的,与组匹配且相互兼容的片段之间采用连通的图形式表示(overlapgraph),每个片段在连通图内采用一个节点（Node）表示,连通图的边线位于每对兼容的片段之间,黄色，蓝色与红色片段是来源于独立的异构体，但其他片段可能来源于3个异构体的任意一个,异构体（Isoform）则可利用的片段拼接起来,其次，Cufflinks软件采用DlorthsTheorem理论重构可能出现的异构体模型。最后，转录本的丰度，即表达量，被计算和评估,Cufflinks软件是采用一种统计学模型计算观测到的每个片段的概率是与其来源的转录本的表达量成线性函数关系,因为每段序列的双端是被测图3.1.1Cufflinks组装原理图示表3.1.1 GTF格式说明第4 第5 第6列 得分和数字，是注释信息可能性的说明,”.”表示为空第7 序列的方向，+代表正链，-代表负链，.代表未第9 注释信息，包括编号，转录本编号，表达量针对所有样本组装且再次组装的Cufflinks组装转录本，采用 pare与已知pare主要给出了5种类型的输出文件，5种类型文件的详细说明见http://cole-trapnell-lab..io pare/#transfrag-class-codes，现简要给表 Potentiallynovelisoform(fragment):atleastonesplicejunctionissharedwitha Atransfragfallingentirelywithinareferenceintron（与已知内含子链特异完全匹配的转录本） Genericexonicoverlapwithareferencetranscript（与已知多个外显子复杂的转录本）Singleexontransfragoverlapareferenceexonandatleast10bpofareference indicatingapossiblepre-mRNAfragment.（与已知单个外显 sAnsAnintronofthetransfragoverlapsareferenceintronontheoppositestrand(likelydueto transcrip)3.2.1表 tmap文 =AABR.1ENSRNOT00000087543ENSRNOG00000040316ENSRNOT00000062051=AABR.1ENSRNOT00000062051ENSRNOG00000023659ENSRNOT00000028992=Raet1dENSRNOT00000028992ENSRNOG00000042852ENSRNOT00000062030=LOC679782 000000000000000000000000000000 Thegene_nameattributeofthereferenceGTFrecordforthistranscript,present.Otherwisegene_idis Thetranscript_idattributeofthereferenceGTFrecordforthistranscript ThetypeofrelationshipbetweentheCufflinkstranscriptsincolumn4andthereferencetranscript(asdescribedintheClassCodessection TheCufflinksinternalgeneid TheCufflinksinternaltranscriptidTheexpressionofthistranscriptexpressedasafractionofthemajorisoform

thegene.Rangesfrom1to TheexpressionofthistranscriptexpressedinFPKM Thelowerlimitofthe95%FPKMconfidenceinterval Theupperlimitofthe95%FPKMconfidenceinterval Theestimatedaveragedepthofreadcoverageacrossthetranscript. Thelengthofthetranscript TheCufflinksIDofthegene’smajorisoform Thelengthofthetranscriptinref1 列说第2 第3 第4 第5 第6列 得分和数字，是注释信息可能性的说明,”.”表示为空第7 序列的方向，+代表正链，-代表负链，.代表未第8 9列 将已知间隔区域内的新转录本利用Blastx程序与Uniprot蛋白质序列数据库分别进行比对，Evalue1E-3。比对结果中Uniprot蛋白质靶标序列的GO注释信息对新转录本进行GO功能分类注释。采用KEGG的KAAS注释工具对新转录本进行KEGGPathway注释。针对新转录本序列的正链和反向互补链，利用Trinity软件中的transdecoder程序基于模型（log似然比，基于可编码与非编码的log似然比）原理预测转录本潜在蛋白intero数据库成员包括Coils、Gene3D、Pfam、PRINTS、ProSitePatterns、ProSiteProfiles、SMART、SUPERFAMILY、TIGRFAM、ProDom、PIR数据库。采用interproscan软件可以对新蛋白质序列通过序列比对或者HMM算法等搜索与interpro蛋白质特征序列匹配预测蛋白质sp|Q05481|ZNF91_HUMANZincfingerprotein91OS=HomosapiensGN=ZNF91PE=2SV=2sp|O43345|ZN208_HUMANZincfingerprotein208OS=HomosapiensGN=ZNF208PE=2SV=2sp|A8MXY4|ZNF99_HUMANZincfingerprotein99OS=HomosapiensGN=ZNF99PE=2SV=3sp|A6NN14|ZN729_HUMANZincfingerprotein729OS=HomosapiensGN=ZNF729PE=2SV=4sp|P17038|ZNF43_HUMANZincfingerprotein43OS=HomosapiensGN=ZNF43PE=2SV=4sp|Q8TF20|ZN721_HUMANZincfingerprotein721OS=HomosapiensGN=ZNF721PE=2SV=2sp|Q8NB50|ZFP62_HUMANZincfingerprotein62homologOS=HomosapiensGN=ZFP62PE=1SV=3sp|Q96IR2|ZN845_HUMANZincfingerprotein845OS=HomosapiensGN=ZNF845PE=2SV=3sp|Q6ZNA1|ZN836_HUMANZincfingerprotein836OS=HomosapiensGN=ZNF836PE=2SV=2sp|Q9HCG1|ZN160_HUMANZincfingerprotein160OS=HomosapiensGN=ZNF160PE=2SV=3 E E 查询序列IDQueryLength SbjctLength QueryAlignment SbjctAlignment EValue cellATPmotoraxonemaldyneinciliumorflagellum-dependentcelldyneinnucleotide 序列ID GeneOntologyID GO名称 GO类别3.3.3 NON-ALCOHOLICFATTYLIVERDISEASE(NAFLD)ALZHEIMER'SDISEASEPARKINSON'SDISEASEComplementandcoagulationcascades NOD-LIKERECEPTORSIGNALINGPATHWAYOthertypesofO-glycanbiosynthesis Metabolism/GlycanbiosynthesisandmetabolismCELLADHESIONMOLECULES(CAMS) 序列ID KEGGPathwayID KEGGPathway类别 KEGGPathway名1 .1 .1 .3 . .1 .1 . .1 . . 第1 第2 第3列 区域类型genemRNAexonCDSfive_primer_UTRthree_primer_UTR，其中CDS代表潜在蛋白质序列区域，*_UTR代表潜在UTR序列区域第4 第5 第7 第9列 编码蛋白序列区域编号，长度，方向(+/-)，类型，其中分为internal：不包括起 子的不完整蛋白质序列； 子的不完整蛋白质序列；complete：包括起 3.3.5interproscanZincfingerC2H2ZincfingerC2H2ZincfingerC2H2ZincfingerC2H2ZincfingerC2H2ZincfingerC2H2Zincfinger,C2H2Zincfinger,C2H2ZincfingerC2H24Zincfinger,Zincfinger,Zincfinger,Zincfinger,Zincfinger,Zincfinger,Zincfinger,Zincfinger,-Zincfinger, 蛋白质序列长度 Interpro数据库的成员

Interproscan预测算法各个成员预测概率值 Interproscan编号 结果输出文件路径在每个样本的./ 所有样本的整理表达水平比较分析进行采用表达量(normalize后画BoxplotPCA。Boxplot图是利用数据中的五个统计量：最小值、第一四分位数（25%）、中位数（50%）、计算(pearsoncorrelation)相关系数R2。样本相关系数聚类图是通过计算样本间相关系数构建的层级聚类图。PCA图是通过降维度的方式来看样本间的相似程度，样本用不同图 已知转录本表达量丰度Boxplot4.1.24.1.3图 样本的PCA针对所有样本采用HTSeq软件统计的原始序列计数，针对有重复样本的实验设计利用DESeq软件对不同样本组之间筛选差异表达的已知，根据客户要求，分别S1*vsS2*；S1*vsS3*；S1*vsS4*；S2*vsS3*；S2*vsS4*;S3*vsS4*六处我们选择满足大于等于Pvalue<=0.05差异表达范围筛选两组之间的差异。下面展示的是S1*vsS2*组的部分结果。表 .表 表 meannormalisedcountsfrom meannormalisedcountsfrom meannormalisedcounts,averagedoverallsamplesfrombothconditions thelogarithm(tobasis2)ofthefoldchange

pvalueforthestatisticalsignificanceofthispvalueadjustedformultipletestingwiththeBenjamini-Hochbergprocedure(seetheRfunctionp.adjust),whichcontrolsfalsediscoveryrate(FDR)图 log2(foldchange)，纵坐标代表-log10(Pvaluefoldchange2差异，绿色横线对应Pvalue=0.05红色是差异区域 图 不同差异集venn图 差异采用相关系数(Pearsoncorrelationcoefficient)且满足R2>=0.9来预测共表达对从而构建共表达网络。由于比较多，以下只显示部分，客户可用cytoscape软件将差异对整个网络文件导入根据需要进行网络可视化。表 图 差异表达GO功能分采用clusterProfiler软件进行GO功能分析，GO功能分析是针对全/转录本和差异/转录本进行功能注释和归类，如果从全/转录本筛选到差异/转录本列表，可以对差异/转录本列表进行GO功能富积分析，GO功能富积分析的方法：将全部/转录本作为背景列表，差异/转录本列表作为从背景列表中筛选出来的候选列表，利用FisherGO功能集在差异/转录本列表中是否显著富积的PP值经Benjamini&Hochberg多重检验纠正。表 fattyacidmetabolicsmallmoleculealcoholmetabolicsterolmetaboliccarboxylicacidorganichydroxymetabolicorganicacidsmallmoleculesecondaryalcoholcholesterolmetabolicGOGO

注释在这个GO的总 q ， 图 显著富积GO柱状图 显著富积GO树状 Pathway Pathway注释在这个pathway的 列表中pathway PBenjamini&HochbergP 个的个的

q候 GOp-差异表达的KEGGpathway功能分KEGGpathway功能分析是针对全/转录本和差异/转录本进行KEGG数据库PathwayKEGGpathwayGO表 KEGGPathway功能富积分析结果文CarbonSteroidFc BcellreceptorsignalingGlyoxylateandBiosynthesisofunsaturatedfattyPPARsignalingPyruvate图 显著富积KEGGpathway柱状图5.2.2差异KEGG富集散点富集程度通过Richfactor、Qvalue和富集到此通的个数来衡量。其中Richfactor指差异表达的中位于该pathway条目的总数。Richfactor越大，表示富集的程度越大。qvalue是做过多重假设检验校正之后的Pvalue，越接近于零，表示富集越显著。图 显著富积KEGGpathway从STRING蛋白质相互作用网络数据库()中查询和该物种的蛋白质相互作用网络，从而构建差异表达蛋白质相互作用网络，如果物种在STRING数据库中不存在，那么近源物种的蛋白质相互作用网络，采用Blastx或者Blastp将该物物种的蛋白质相互作用网络，映射的方法是如果该物种A的近源物种同源列表中的1个或者多个与该物种B的近源物种同源列表中的1个或者多个存在蛋白质相互作用关系，那么认为A与B存在可能的蛋白质相互作用，从而构建差异蛋白质相互作用网络。由于差异数目较多，选取度排序前50个来绘制蛋白质相互作用网络图。客户可用cytoscape软件将差异对整个网络文件导入根据需要进行网络可 图5.3.1差异蛋白质相互作用网络注：strings_network.txt是该物种差异蛋白质相互作用网络文件，第1列是A对应的蛋白，第2列是B对应的蛋白，第3列是A和B的相互作用得分，第4列是A对应的entrezID，第5列是A对应的symbol，第6列是B对应的entrezID，第7列是Bsymbol结果输出文件路径./5_GSEASNP和INDELSNP和INDEL采用BWA软件对预处理序列与物种的参考组序列进行序列比对，针对比对结果文件采用picard-tools去掉由于PCR过程的重复序列，采用GTAK3分别进行SNP和INDEL表 每个样本的SNP或者INDEL分析结果VCF格式文1.TC1.TG1.TA1.GCS1.GAS1.TCS1 A1 AS1 GS1 AS1 AS variant FILTER分数与覆盖深度比值过滤，FSFilter表示参考碱基和变异碱基的序列链方向分布偏见过滤，ReadPosRankSumFilter3置分布偏见过滤，HaplotypeScoreFilter表示基于单倍型得分过滤

变异碱基相关信息,过滤指标的数值.ABHom/ABHet代表变异的纯合和杂合；ABHom:AlleleBalanceforhets(ref/(ref+alt));ABHom:AlleleBalanceforhoms(A/(A+O));AC(AlleleCount)表示该Allele的数目；AF(AlleleFrequency)表示AlleleAN(AlleleNumber)表示AlleleBaseQRankSum:Z-scorefromWilcoxonranksumtestofAltVs.Refbasequalities；DB:覆盖碱基的深度;Dels:FractionofReadsContainingSpanningDeletions;FS:Phred-scaledp-valueusingFisher'sexacttesttodetectstrandbias;HaplotypeScore:Consistencyofthesitewithatmosttwosegregatinghaplotypes;MLEAC:umlikelihoodexpectation(MLE)fortheallelecounts;MLEAF:umlikelihoodexpectation(MLE)fortheallelefrequency;MQ：RMSMapQuality；MQ0：TotalMapQualityZeroReads；MQRankSumZ-scoreFromWilcoxonranksumtestofAltRefreadmapqualities；OND:Overallnon-diploidSample型型）:AD（REFALT深(alleles/(alleles+non-alleles));Sample型型）:AD（REFALT深度）:GQ（代表型准确度的质量分数）:PL（0/0纯合,0/1杂合突变,1/1纯合突变这三种型的值用来判断GQ，值越大错误越大）FORMAT对应的数值。如0/1:20,4:24:97:97,0,418：0/1表示杂合,其中REF20个，ALT为4个，深度为24个，型的准确性97，（0/0,0/1,and1/1的值为97,0,418由于0为错误最小值所以表型为0/1）SNP和INDEL采用ANNOVAR对所有样本的SNP和INDEL位点进行 水平(gene-based),区域水平(region-based)和已知位点水平(filter-based)水平的注释。 (ncRNA_exonic)、非编码蛋白质的外显子区域(ncRNA_intronic)、间隔区域(intergenic)、可变性剪切位点(splicing)、UTR5区域(UTR5)、UTR3区域(UTR3)、上游区域(upstream)、下游区域(downstream)、编码蛋白质的外显子区域内是否改变氨基酸(nonsynonymous、synonymous、stopgain、stoploss、frameshiftinsertion、frameshiftdeletion)编码；区域特征水平的注释包括跨7个物种的保守性区域(phastConsElements7way)、组的大片段重复区域(genomicSuperDups)、染色体遗传标记(cytoBand)、CpG岛、DGV(DatabaseofGenomicVariants)GWASfilter水平的注释包括1000GenomesProjectannotationsdbSNPannotations、ESP(exomesequencingproject)annotations、CG(completegenomics)frequencyannotationspopulationfrequencyensembleannotationsCLINVARannotations、COSMICannotations。表 每个样本的SNP或者INDEL注释结果文件（其他信息见结果文件1TC1TG1TA1GC1GA1TC1TA1AG1CA1CA

定位到exonic,intergenic,intronic,UTR5,UTR3,splicing,ncRNA_intronicncRNA_splicingncRNA_UTR3ncRNA_UTR5GeneDetai ExonicFun

输出格式为 ID:转录本ID:外显子编号:变异位点，例如OR4F5: :exon1:c.A421G:p.T141A，转录本 表 1 1 1 1 1 1 1 1 1 1 CC.................................................................. variant型，A/G表示A是ref碱基，G是alt碱基，.表示型缺结果输出文件路径 类型包括SKIP、XSKIP、XMSKIP、IR、XIR、MIR、XMIR、AE和XAE。表表

Alternative5frststartAlternative3lasexon(transcriptterminalsite)

单个外显子跳跃，ON代表包含，OFF SkippedApproximate

单个外显子跳跃，ON代表包含，

Multi-exonSKIP(MSKIP_ON,

多个外显子跳跃，ON代表包含， Intronretention(IR_ON, 单个内含子延伸，ON代表延伸， ApproximateIR(XIR_on, 单个内含子延伸，ON代表延伸， Multi-IR(MIR_ON, Alternativeexonend(5,3,

多个内含子延伸，ON代表延伸，代表不延伸，内含子边界模糊53 Approximate 53 111111 -0-0-0-0+0+0++

可变性剪切特征(forTSS,TTS-insideboundaryofalternativemarginalexon;for*SKIP_ON,thecoordinatesoftheskippedexon(s);forcoordinatesofallthe