在这一 部分中,对前面生成的SAM文件进行定量,用到的工具是htseq
,如下所示:
|
|
生成的数据如下所示:
用Notepad++打开其中的一个文件,例如SRR3589959count
,如下所示:
数据分为两列,第1列是基因名,第2列是reads数。现在就要把这4个文件给合并起来,构成一个矩阵,如下所示
|
|
在这一 部分中,对前面生成的SAM文件进行定量,用到的工具是htseq
,如下所示:
|
|
生成的数据如下所示:
用Notepad++打开其中的一个文件,例如SRR3589959count
,如下所示:
数据分为两列,第1列是基因名,第2列是reads数。现在就要把这4个文件给合并起来,构成一个矩阵,如下所示
|
|
- 搞懂hisat2的用法。
- 直接去hisat2的主页下载index文件即可,然后把fastq格式的reads比对上去得到sam文件。
- 接着用samtools把它转为bam文件,并且排序(注意N和P两种排序区别)索引好,载入IGV,再截图几个基因看看。
- 顺便对bam文件进行简单QC,参考直播我的基因组系列。 来源于生信技能树:http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
HISAT2:http://ccb.jhu.edu/software/hisat2/index.shtml
参考资料:http://blog.biochen.com/archives/337
STAR:https://codeload.github.com/alexdobin/STAR/zip/master
参考资料:http://www.bio-info-trainee.com/727.html
TopHat:http://ccb.jhu.edu/software/tophat/index.shtml
参考资料:http://blog.sina.com.cn/s/blog_8808cae20101amqp.html
RapMap:https://github.com/COMBINE-lab/RapMap
参考:https://academic.oup.com/bioinformatics/article/32/12/i192/2288985/RapMap-a-rapid-sensitive-and-accurate-tool-for
CIDANE:http://ccb.jhu.edu/software/cidane/
参考文献:https://genomebiology.biomedcentral.com/articles/10.1186/s13059-015-0865-0
CLASS2 :https://sourceforge.net/projects/splicebox/files/?source=navbar
参考文献:https://academic.oup.com/nar/article/44/10/e98/2516329/CLASS2-accurate-and-efficient-splice-variant
内容主要参考:http://www.360doc.com/content/16/1223/13/29483982_617058719.shtml
人类和小鼠的索引有现成的,HISAT2官网可以直接下载进行序列比对。如下图所示:选择hg19和mm10的index,文章中RNA-Seq测序数据,可以包括人类的3个数据和小鼠的4个数据,因此需要小鼠和人类的索引。
|
|
解压后的结果如下所示(小鼠基因组的index):
需要注意的是,小鼠的基因组index这几个件,前面都是genome
为前缀。
Usage:hisat2 [options]* -x {-1 -2 | -U | —sra-acc } [-S ]
参数:
-x 指定index文件
-1 双端测序第一个文件
-2 双端测序第二个文件
-U 单端测序文件
—sra-acc SRA accession number
-S 指定输出的格式,一般指定为sam
代码如下:
这里需要注意的是,在进行比对的时候,我们使用的命令是这个样子的/media/w/新加卷/Download/reference/mm10/genome
,最后一个是genome
,它代表的是前面的小鼠基因组index的8个文件,在输入的时候,不能只输入文件夹,也就是只输入/media/w/新加卷/Download/reference/mm10
这一部分,还要添加上genome
。
序列比对后得到SAM文件,这里只看小鼠数据的比对结果(也就是59到62这一部分),如下所示:
SAM(sequence Alignment/mapping)数据格式是目前高通量测序中存放比对数据的标准格式,当然他可以用于存放未比对的数据。所以,SAM的格式说明
而目前处理SAM格式的工具主要是SAMTools,这是Heng Li大神写的,除了C语言版本,还有Java的Picard,Python的Pysam,Common lisp的cl-sam等其他版本。SAMTools的主要功能如下:
view: BAM-SAM/SAM-BAM 转换和提取部分比对
sort: 比对排序
merge: 聚合多个排序比对
index: 索引排序比对
faidx: 建立FASTA索引,提取部分序列
tview: 文本格式查看序列
pileup: 产生基于位置的结果和 consensus/indel calling
最常用的三板斧就是格式转换,排序,索引。而进阶教程就是看文档提高。
|
|
最终的BAM文件结果如下所示:
在UCSC下载hg19参考基因组,群主博客有详细说明,从gencode数据库下载基因注释文件,并且用IGV去查看你感兴趣的基因的结构,例如TP53,EGFR等等。 截图几个基因的IGV可视化结构!还可以下载ENSEMBL,NCBI的GTF,也导入IGV看看,截图基因结构。了解IGV常识。 来源于生信技能树:http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
测序得到的是几百bp的短read, 相当于把拼图打散了给你。如果没有参考基因组,从头(de novo)组装等于是重走人类基因组计划的老路,也就是打散了拼图,却不告诉你原来是什么样子,那么任务将会及其艰巨。还好人类基因组已经组装好了,我们只需要把我们测得序列回贴(mapping)回去,毕竟人与人之间的差距只有不到1%差异, 允许mismatch就行。
参考基因组对新手来说,是一个很大的坑,hg19、GRCH37、 ensembl 75这3种基因组版本应该是大家见得比较多的了,国际通用的人类参考基因组,其实他们储存的是同样的fasta序列,只是分别对应着三种国际生物信息学数据库资源收集存储单位,即NCBI,UCSC及ENSEMBL各自发布的基因组信息而已。有一些参考基因组比较小众,存储的序列也不一样,比如BGI做的炎黄基因组,还有DNA双螺旋结构提出者沃森(Watson)的基因组,还有2016年发表在nature上面的号称最完善的韩国人做的基因组。前期我们先不考虑这些小众基因组,主要就下载hg19和hg38,都是UCSC提供的,虽然hg38相比hg19来说,做了很多改进,优点也不少,但因为目前为止很多注释信息都是针对于hg19的坐标系统来的,我们就都下载了,正好自己探究一下。也顺便下载一个小鼠的最新版参考基因组吧,反正比对也就是睡个觉的功夫,顺便分析一下结果,看看比对率是不是很低。
这里下载的是USCS版本的参考基因组。
|
|
由于下载速度慢,我就改用迅雷下载了,将所有染色体信息整合后,得到新的hg19.fa文件。
参考基因组是一部无字天书,要想解读书中的内容,需要额外的注释信息协助。因此下载完参考基因后,就是去gencode数据库(http://www.gencodegenes.org/)下载基因组注释文件。简单来讲注释文件就是基因组的说明书,告诉我们哪些序列是编码蛋白的基因,哪些是非编码基因,外显子、内含子、UTR等的位置等等。注释文件在以上三个提供参考基因组的网站中都有提供,比如Ensemble。但是现在最权威的人类和小鼠基因组的注释还属Gencode数据库。
进入人和小鼠基因组注释信息官网GENCODE,选择data->human->GRCh37-mapped Releases,下载最新第26版本的hg19人类基因组注释信息。点击进入下载页面,将GTF和GFF3全部下载,解压,如下所示:
可以看到已经有新的参考基因组注释了,不过现在还用老的,跟参考资料中的保持一致一样,单击26,进入下载界面,下载第一行的注释信息,即Comprehensive gene annotaton中的GTF和GFF3,如下所示:
GFF全称为general feature format,这种格式主要是用来注释基因组。
GTF全称为gene transfer format,主要是用来对基因进行注释。
GFF文件是一种用来描述基因组特征的文件,现在我们所使用的大部分都是第三版)(GFF3)。GFF允许使用#作为注释符号,例如很多GFF文件都会使用如下的两行来表明其版本其创建日期:
GFF文件每一列所代表的含义后面表格中有,但请注意,它的第3列feature type是不受约束的,你可以使用任意的名称。我们需要注意的是GFF文件的第9列,从第二版开始(GFF2),所有的属性都以标签=值的方式呈现,各个属性之间以;作为分隔符
在最新版本的GFF文件中(GFF3),有一些是已经预先定义的属性特征,并且这些特征往往还有特殊的含义:ID这个标签实在各行都要有的;另外有一个Parent的属性,它表明了当前的特征是Parent特征的子集。
|
|
当前所广泛使用的GTF格式为第二版(GTF2),它主要是用来描述基因的注释。GTF格式有两个硬性标准:
列|GTF2|GFF3
|:-:|:-:|:-:|
|1.reference sequence name|same|same|
|2.annotation source |same| same|
|3.feature type| feature requirements depend on software| can be anything|
|4.start coordinate| same| same|
|5. end coordinate |same |same|
|6.score |not used |optional
|7.strand |same| same|
|8.frame |same |same|
|9.attributes |空格分隔|=分隔|
Integrative Genomics Viewer(IGV)是一种探索大型综合基因组数据的高性能交互式可视化工具。它支持各种各样的数据类型,包括基于芯片测序、二代测序数据和基因组注释数据等。
|
|
如下所示:
需要用安装好的sratoolkit把sra文件转换为fastq格式的测序文件,并且用fastqc软件测试测序文件的质量! 作业,理解测序reads,GC含量,质量值,接头,index,fastqc的全部报告,搜索中文教程,并发在论坛上面。 来源于生信技能树:http://www.biotrainee.com/forum.php?mod=viewthread&tid=1750#lastpost
下载的原始数据是sra格式,SRA(Sequence ReadArchive)数据库是用于存储二代测序的原始数据,包括 454,Illumina,SOLiD,IonTorrent,Helicos 和 CompleteGenomics。除了原始序列数据外,SRA现在也存在raw reads在参考基因的比对信息。
Fastq格式是一种基于文本的存储生物序列和对应碱基(或氨基酸)质量的文件格式。最初由桑格研究所(Wellcome Trust Sanger Institute)开发出来,现已成为存储高通量测序数据的事实标准。
FASTQ文件中每个序列通常有四行:
|
|
如何判断是Phred64 还是 Phred33 ? ASCII值小于等于58(相应的质量得分小于等于25)对应的字符只有在Phred+33的编码中被使用,所有Phred+64所使用的字符的ASCII值都大于等于59。在通常情况下,ASCII值大于等于74的字符只出现在Phred+64中。如果是最近两年的测序数据,一般都是Phred33形式的。参考文章:http://blog.csdn.net/huyongfeijoe/article/details/51613827
将sra格式转换为fastq格式的工具是sratoolkit的fastq-dump,可以先用fastq-dump -h看一下帮助文件,分为如下几个部分:
输入: -A|—accession 序列号
处理中: Read Splitting, Full Spot Filters, Common Filters, Filters based on alignments, Filters for individual reads。 基本都是些过滤参数。不太常用
输出: -O|—outdir 输出文件夹, -Z|—stdout 输出到标准输出, —gzip/—bzip2 输出为压缩格式
多文件选项: 常用的就是—split-3,split-3表示如果是单端测序则一个sra文件出来一个fastq文件,如果是双末端,则一个sra问件对应两个fastq问件SRRXXXXXX_1.fastq,SRRXXXXXX_2.fastq
格式化: 分为序列, 质量等,不常用
—gzip 使得输出的结果是.gz 的格式
现在将disk2/sra下的sra文件转化为fastq格式,其中输入目录为disk2/sra,输出目录为/disk2/data/rna-seq,命令如下所示:
|
|
转换后的结果如下所示:
用法如下所示:
fastqc [-o output dir] [—(no)extract] [-f fastq|bam|sam] [-c contaminant file] seqfile1 .. seqfileN
- 参数: -o 输出目录,需自己创建目录
- —(no)extract 是否解压输出文件,默认是自动解压缩zip文件。加上—noextract不解压文件。
- -f 指定输入文件的类型,支持fastq|bam|sam三种格式的文件,默认自动识别。
- -t 同时处理的文件数目。
- -c 是contaminant 文件,会从中搜索overpresent 序列。
代码如下所示:
|
|
质控结果有14个html文件,可以选择用浏览器打开查看最终的QC reports。随便打开一个质控结果文件,如下所示:
左边是目录概要,一共12项内容,如下所示:
可以点击想要看的结果,右边会跳转到特定详细的可视化结果。绿色代表“通过”,黄色代表“警告”,红色代表“不通过,失败”,如下所示:
基本的数据统计包括文件名,文件类型,编码形式,总的序列数,质量差的序列,序列平均长度,GC含量。
图片太大,下部图片如下所示:
上图表示的是每个read各位置碱基的测序质量,具体的理解如下所示:
一般要求此图中,所有位置的10%分位数大于20,也就是我们常说的Q20过滤(即质量数低于20的reads清除掉)W。
检查reads中每一个碱基位置在不同的测序小孔之间的偏离度,蓝色代表偏离度小,质量好,越红代表偏离度越大,质量越差。纵轴是tail的Index编号。
这个图主要是为了防止,在测序过程中,某些tail受到不可控因素的影响而出现测序质量偏低。蓝色代表测序质量很高,暖色代表测序质量不高,如果某些tail出现暖色,可以在后续分析中把该tail测序的结果全部都去除
质量的分布,当峰值小于27时,警告;当峰值小于20时,fail。可以看出这个测序的报告峰值在38左右。
假如我测的1条序列长度为101bp,那么这101个位置每个位置Q之的平均值就是这条reads的质量值。该图横轴是0-40,表示Q值,纵轴是每个值对应的reads数目
对所有reads的每一个位置,统计ATCG四种碱基的分布,横轴为位置,这里是1-51bp,纵轴为碱基含量(百分数),正常情况下每个位置每种碱基出现的概率是相近的,四条线应该平行且相近。当部分位置碱基的比例出现bias时,即四条线在某些位置纷乱交织,往往提示我们有overrepresented sequence的污染。本结果前10个位置,每种碱基频率有明显的差别,说明有污染。当任一位置的A/T比例与G/C比例相差超过10%,报”WARN”;当任一位置的A/T比例与G/C比例相差超过20%,报”FAIL”。
理论上来说,A和T应该相等,G和C应该相等,但是一般测序的时候,刚开始测序仪状态不稳定,很可能出现上图的情况。像这种情况,即使测序的得分很高,也需要cut开始部分的序列信息,一般遇到到这种情况,会cut前面5bp。
统计reads的平均GC含量的分布。横轴是0 - 100%; 纵轴是每条序列GC含量对应的数量。红线是实际情况,蓝线是理论分布(正态分布,均值不一定在50%,而是由平均GC含量推断的),两个应该比较接近才比较好。 曲线形状的偏差往往是由于文库的污染或是部分reads构成的子集有偏差(overrepresented reads)。形状接近正态但偏离理论分布的情况提示我们可能有系统偏差。偏离理论分布的reads超过15%时,报”WARN”;偏离理论分布的reads超过30%时,报”FAIL”。当红色的线出现双峰,基本肯定是混入了其他物种的DNA序列。这张图中的信息良好。
当测序仪器不能辨别某条reads的某个位置到底是什么碱基时,就会产生“N”,统计N的比率。正常情况下,N值非常小。当任意位置的N的比例超过5%,报”WARN”;当任意位置的N的比例超过20%,报”FAIL”。
reads长度分布,每次测序仪测出来的长度在理论上应该是完全相等的,但是总会有一些偏差。比如此图中,51bp是主要的,但是还是有少量的50和52的长度,不过数量比较少,不影响后续分析。当测序的长度不同时,如果很严重,则表明测序仪在此次测序过程中产生的数据不可信。当reads长度不一致时报”WARN”;当有长度为0的read时报“FAIL”。
统计不同拷贝数的reads的频率。测序深度越高,越容易产生一定程度的duplication,这是正常的现象,但如果duplication的程度很高,就提示我们可能有bias的存在。横坐标是duplication的次数,纵坐标是duplicated reads的数目,以unique reads的总数作为100%。下图中,横坐标表示序列重复水平;纵坐标代表重复水平序列占所有序列的百分比。下图中大于10个重复的reads占总序列的20%以上,其他依次类推。当非unique的reads占总数的比例大于20%时,报”WARN”;当非unique的reads占总数的比例大于50%时,报”FAIL“。
一条序列的重复数,因为一个转录组中有非常多的转录本,一条序列再怎么多也不太会占整个转录组的一小部分(比如1%),如果出现这种情况,不是这种转录本巨量表达,就是样品被污染。这个模块列出来大于全部转录组1%的reads序列,但是因为用的是前200,000条,所以其实参考意义不大,完全可以忽略。
接头含量。①此图衡量的是序列中两端adapter的情况;②如果在当时fastqc分析的时候-a选项没有内容,则默认使用图例中的四种通用adapter序列进行统计;③本例中adapter都已经去除,如果有adapter序列没有去除干净的情况,在后续分析的时候需要先使用cutadapt软件进行去接头。
①这个图统计的是,在序列中某些特征的短序列重复出现的次数;② 我们可以看到1-8bp的时候图例中的几种短序列都出现了非常多的次数,一般来说,出现这种情况,要么是adapter没有去除干净,而又没有使用-a参数;要么就是序列本身可能重复度比较高,如建库PCR的时候出现了bias;④ 对于这种情况,我的办法是可以cut掉前面的一些长度,可以试着cut 5~8bp
|
|
本系列课程学习的文章是:AKAP95 regulates splicing through scaffolding RNAs and RNA processing factors. Nat Commun 2016 Nov 8;7:13347. PMID: 27824034 很容易在文章里面找到数据地址GSE81916 这样就可以下载sra。
人类,小鼠;
高通量测序
pre-mRNA的可变剪切与高等生物基因表达的复杂性有关,但是可变剪切位点的选择却未研究清楚。作者在以前的研究中发现了一个染色质相关的蛋白AKAP95(AKAP8)在强化染色质转录方向有着明显的影响。在本文的研究中,APAKP95能与许多涉及转录和RNA加工的因子相互作用,调节pre-mRNA的剪接。AKAP95直接促进内源基因FAM126A的特异性外显子的体外剪接和包含。 AKAP95N-末端富含YG的结构域在与RNA加工因子(例如hnRNP蛋白)结合的过程中有着重要的作用,它的锌指结构域对于pre-mRNA的结合至关重要。基因组结合的分析显示,AKAP95优先结合人类转录组中大量pre-mRNA的近端内含子区域,AKAP95的敲除能导致多外显子增加的减少。AKAP95还可以选择性地协调hnRNP H/ F和U蛋白来调节可变剪接事件。进一步表明,AKAP95与自身直接相互作用。总之,实验的研究结果表明,AKAP95是pre-mRNA剪接的新型,并且是正向调节因子,它有可能是一个转录和剪接调节整合子,并且实验提出了一个模型,即AKAP95通过支架RNA和RNA加工因子调节pre-mRNA剪接,并促进剪接位点之间的信息交流。
原英文摘要:
Alternative splicing of pre-mRNAs significantly contributes to the complexity of gene expression in higher organisms, but the regulation of the splice site selection remains incompletely understood. We have previously demonstrated that a chromatin-associated protein, AKAP95 (AKAP8), has a remarkable activity in enhancing chromatin transcription. In this study, we have shown that AKAP95 physically interacts with many factors involved in transcription and RNA processing, and functionally regulates pre-mRNA splicing. AKAP95 directly promotes splicing in vitro and the inclusion of a specific exon of an endogenous gene FAM126A. The N-terminal YG-rich domain of AKAP95 is important for its binding to RNA processing factors including selective groups of hnRNP proteins, and its zinc finger domains are critical for pre-mRNA binding. Genome-wide binding assays revealed that AKAP95 bound preferentially to proximal intronic regions on a large number of pre-mRNAs in human transcriptome, and AKAP95 depletion predominantly resulted in reduced inclusion of many exons. AKAP95 also selectively coordinates with hnRNP H/F and U proteins in regulating alternative splicing events. We have further shown that AKAP95 directly interacts with itself. Taken together, our results establish AKAP95 as a novel and mostly positive regulator of pre-mRNA splicing and a possible integrator of transcription and splicing regulation, and support a model that AKAP95 regulates pre-mRNA splicing via through scaffolding RNAs and RNA processing factors and facilitating the splice site communication.
样本1-8是RIP-seq的数据,研究的是AKAP95与转录组的结合,样本9-15是mRNA-seq的数据,研究的是在293细胞或小鼠ES细胞中敲减了AKAP95的数据(其中9-11是293细胞的数据,有3个,12-15是小鼠ES细胞的数据,有4个)。
1.AKAP95的N末端区域参与该蛋白与RNA加工因子的结合,这些RNA加工因子包括hnRNPs。
2.APAK95与许多RNA加工因子的结合暗示了AKAP95在RNA加工方面发挥了作用,这些作用包括pre-RNA剪接。实验结果表明,AKAP95直接调控minigene pre-mRNA剪接,AKAP95的N末端与ZFs是其功能发挥的结构域。
3.为了寻找AKAP95剪接调控的内源靶点。经过一系列实验表明:AKAP95和hnRNP F参与调节FAM126A外显子11的剪接(实验中提到了敲减AKAP95后,FAM126A的外显子11会被跳过,而在野生型中则正常剪接)。
4.利用深度测序的方法寻的找内源AKAP95与人类转录组的结合。测序手段是RIP,其结果是发现AKAp95能与pre-mRNAs结合。
5. (核心,也就是这个笔记学习的部分)问题:评估AKAP95对可变剪接的全局影响,实验手段是:①人类294细胞和小鼠ES细胞中敲减了AKAP95基因->②RNA-seq->③DEXseq分析。AKAP95的敲减会导致更多的外显子使用(exon)减少,这表明AKAP95会全局促进外显子增加(exon inclusion)。看文献我的理解是,有了AKAP95基因后,该基因的表达产物会导致FAM126A的外显子10,11,12连接起来,如果敲低了AKAP95这个基因,FAM126A的外显子10与12连接起来,会跳过外显子11。因此AKAP95的功能就在于调节FAM126A的可变剪接。exon inclusion我的理解是,外显子的利用,以此文献为例说明,exon inclusion指的就是FAM126A的外显子11的利用,AKAP95会促进AKAP95外显子11的利用,如果没有AKAP95,则FAM126A的外显子10就会跳过外显子11,与外显子12连接起来,exon inclusion就降低。
6.后续的实验作者又研究了AKA95的一些作用机理,略。
a图是火山图。火山图常用于展示显著差异表达的基因,这里有两个关键词:显著是指P<0.05,差异表达一般我们按照fold change(倍数变化)="">=2.0作为标准。
当我们拿到基因表达的P值和倍数后,为了用火山图展示结果,一般需要把倍数进行Log2的转化(下图a的x轴),比如某基因在实验组表达水平是对照组的4倍,log2(4)=2,同样的如果是1/4,也就是0.25,转换后的结果就是-2。
P值进行-log10的转化,-log10(0.01)约等于2,由于P值越小表示越显著,所以我们进行-log10(P value)转化后,转化值越大表示差异约显著,比如-log10(0.001)=3 > -log10(0.01)=2 > -log10(0.05)=1.30(下图a的y轴)。
火山图有三种颜色,黄色与红色都是Padj小于0.01的基因(P值经过了FDR校正,即false discovery rate)。红色又表示基因变化超过2倍的基因。b图是饼图,它表示外显子使用(exon usage)的增加与降低数目。0.05,差异表达一般我们按照fold>
在后文的Bioinformatic analyses部分
比对软件:TopHat (v2.0.13)
参考基因组:human reference genome (GRCh37/hg19)
GTF文件: GTF version GRCh37.70
只保留MQ >30的map结果
Picard-tools (v1.126): 计算平均插入大小(mean insert sizes)和标准差
read count: 软件:HTSeq v0.6.0
DESeq (v3.0)
DEXSeq (v3.1)
GO分析:DAVID (http://david.ncifcrf.gov/)
进入NCBI的GEO数据库https://www.ncbi.nlm.nih.gov/geo/
搜索GSE81916,如下所示:
输入GSE81916,单击Search即可,见到数据说明及下载页面,如下所示:
下面有一段话是对这些数据的描述,即:
>
Overall design Samples 1-8 are RNA-immunoprecipitation > (RIP)-seq to determine AKAP95 binding to the
transcriptome. Samples 9-15 are mRNA-seq to determine
effect of AKAP95 knockdown in human 293 cells (9-11) or > mouse ES cells (12-15).
大意就是样本9-15是mRNA-seq的数据,而样本1-8是RNA免疫组化的数据。我们只需要9-15的数据,单击红框链接即可进入数据下载页面,如下所示:
单击第一行的SRP/SRP075/SRP075747右侧的ftp,即可进入数据下载界面:
现在下载SRR3589956-SRR3589962的7个数据(3个人类,4个小鼠),命令如下:
此命令是将数据下载到disk2/sra的目录下,
|
|
当时觉得数据下载慢,我就直接用迅雷下载了,数据下载后如下所示:
转录组的这几篇入门笔记主要是参照了生信技能树论坛,公众号,生信媛公众号,徐更洲的博客,沈梦圆的博客及公众号,Jimmy等人的博客完成的。
先了解一下conda,Anaconda,Miniconda。
conda是一个工具,也是一个可执行命令,其核心功能是包管理与环境管理。包管理与pip的使用类似,环境管理则允许用户方便地安装不同版本的python并可以快速切换。
Anaconda是一个用于科学计算的Python发行版,支持 Linux, Mac, Windows系统,提供了包管理与环境管理的功能,可以很方便地解决多版本python并存、切换以及各种第三方包安装问题。Anaconda则是一个打包的集合,里面预装好了conda、某个版本的python、众多packages、科学计算工具等等,所以也称为Python的一种发行版。
Miniconda则是Anaconda的微缩版,它只包含最基本的内容,即python与conda,以及相关的必须依赖项,对于空间要求严格的用户,Miniconda是一种选择。
bioconda是conda上一个分发生物信息软件的频道。多数软件都可以通过conda来安装,可查看可用软件列表。
|
|
|
|
安装过程基本上就是按Enter或输入yes,仔细看说明就行。
|
|
Conda默认的源访问速度有些慢,可以增加国内的源(国内的源是清华源,其中参考文献中给出的;另外还可以增加几个源,以便于安装更多的软件,尤其是bioconda安装生信类工具。conda-forge通道是Conda社区维护的包含很多不在默认通道里面的通用型软件。r通道是向后兼容性通道,尤其是使用R3.3.1版本时会用到。后添加的通道优先级更高,因此一般用下面列出的顺序添加。
|
|
如果命令行添加不了,可以用vim ~/.condarc
来进行安装,如果安装错误了,可以删除.condarc文件,然后再用命令行来安装。用conda config --get channels
可以查看已经添加的频道,而用conda config --show
可查看已有的配置。如下所示:
|
|
如下所示:
命令是conda search <package ambigious name>
,以搜索fastqc为例:
如下所示:
命令为conda install <package name>
以安装 numpy=1.7.2为例说明
|
|
conda update <package name>
conda remove <package name>
conda -h # 查看conda可用的命令
conda install -h #查看install子命令的帮助
condas可以创建多个分析环境,这是它的优势之一,不过目前用不到,以后用到了再学。
虚拟环境管理:conda比较好用的就是它能够建立多个互不干扰的分析环境。通过conda info --envs
可以查看环境,如下所示:
|
|
目前就只有一个默认环境,也就是root:
|
|
从环境变量中去掉miniconda:打开~/.bash_profile文件,删掉其中miniconda的路径,关闭并保存
删除隐藏的.condarc 、.conda以及.continuum文件
目前所用的生信工具如下所示:
sratoolkit:把NCBI SRA(Sequence Read Archive)数据库中的NGS序列数据从 sra 格式转换到 fastq 格式。
conda安装:
|
|
如果安装不成功,可采用下列安装方式(我的笔记本上上基本上都是conda安装,对于下面的安装方式并没有试过):
|
|
fastqc:二代测序数据质量分析软件。
安装如下:
或者是采用常规的安装方式:
hisats:将测序结果比对到人类参考基因组上。HISAT2是TopHat2/Bowti2的继任者,使用改进的BWT算法,实现了更快的速度和更少的资源占用。
|
|
samtools:处理SAM、BAM文件的工具软件。BAM格式文件是存放高通量测序中比对结果的标准格式文件。
|
|
或者如下安装:
htseq-count:用于reads计数的软件,他能对位于基因组上的一些单位的reads数进行统计,这里所说的单位主要是指染色体上的一组位置区间(我们常见的就是gene exon)
优先conda安装
|
|
如果conda无法安装,采用下列安装方式:
|
|
或者如下:
R:统计与数据可视化工具。
|
|
后面涉及有关R的分析是在Sindows平台下完成的,并没有用到Linux下的R工具。
|
|
在医学研究中有一些现象是难以直接观测的,通常称为不可测现象,它们只能通过其他多个可观测的指标来间接地反映。例如,脑部疾病患者的意识清醒状态是一个不可测现象,但可以通过患者的语言能力、辨识能力、记忆能力、理解能力、思维的逻辑能力等一系列可观测的指标来反映。这里,由于各个可观测指标都不同程度地反映了意识清醒状态这一不可测现象,因此这些可观测指标之间呈现出一定的相关性。很自然地,人们可能认为这些可观测指标之间的相关性主要是由它们反共同反映的不可测现象支配的。
一般地,对于多指标数据中呈现出来的相关性,是否存在对这种相关性起支配作用的潜在因素?如果存在,如何找出这些潜在因素?这些潜在因素是怎样对原始指标起支配作用的?这些问题,都可以通过因子分析来解决。事实上,因子分析是一种从分析多个原始指标的相关关系入手,找到支配这种相关关系的有限个不可观测的潜在变量,并用这些潜在变量来解释原始指标之间的相关性或协方差关系的多元统计分析方法。
案例如下:
例22-2 某医院为了合理地评价该院各月的医疗工作质量,搜集了三年有关门诊人次、出院人数、病床利用率、病床周转次数、平均住院天数、治愈好转率、病死率、诊断符合率、抢救成功率等9个指标数据,如表22-8。现采用因子分析方法,探讨其综合评价指标体系。
|
|
2._config.yml
的配置
添加以下信息:
|
|
参考资料:
主成分分析(principal components analysis,PCA),是由1901年被Pearson首先引入的,1933年由Hotelling作了进一步的发展,主成分分析是从多个数值变量(指标)之间的相互关系入手,利用降维的思想,将多个变量(指标)化为少数几个互不相关的综合变量(指标)的统计方法。
在医学研究中,为了客观、全面地分析问题,常要记录多个观察指标并考虑众多的影响因素,这样的数据虽然可以提供丰富的信息,但同时也使用数据的分析工作更趋复杂化。例如在儿童生长发育的评价中,收集到的数据包括每一个儿童的身高、体重、胸围、坐高、肺活量等十多个指标。怎样利用这类多指标的数据对每一个儿童的生长发育水平作出正常的评价?如果仅用其中任一指标来作评价,其结论显然是片面的,而且不能充分利用已有的数据信息。如果分别利用每一指标进行评价,然后再综合各指标评价的结论,这样做一是可能会出现各指标评价的结论不一致,甚至相互冲突,从而给最后的综合评价带来困难;二是工作量明显增大,不利于进一步的统计分析。事实上,在实际工作中,所涉及的众多指标之间经常是有相互联系和影响的,从这一点出发,希望通过对原始相互关系的研究,找出少数几个综合指标,这些综合指标是原始指标的纯属绵且,它既保留了原始指标的主要信息,且又互不相关。这样一种从众多原始指标之间相互关系入手,寻找少数综合指标以根据原始指标信息的多元统计方法称为主成分分析。
设有m个指标,即$X_{1},X_{2},\dots,X_{m}$,欲寻找可以概括这m个指标主要信息的综合指标$Z_{1},Z+{2},\dots,Z_{m}$,从数学上讲,就是寻找一组常数,即$a_{i1},a_{i2},\dots,a_{im}(i=1,2,\dots,m)$,使这m个指标的线性组合为(以下为公式1-1——:
能够根据m个原始指标$X_{1},X_{2},\dots,X_{m}$的主要信息,且各$Z_{i}(i=1,2,\dots,m)$互不相关,为叙述方便,可以引入以下的矩阵形式,令
则公式1-1可以表示为:
或者是
如果$Z_{1}=\bf a’_{1}X$满足$\bf a’_{1}a_{1}=1$,且
$Var(Z_{1})=\mathop{Max}\limits_{\bf a’a=1}^{}{Var(\bf a’X)}$,刚称$Z_{1}$是原始指标$X_{1},X_{2},\dots,X_{m}$的第一主成分。
通常情况下,如果$Z_{i}=\bf a’X$满足:
(1)$\bf a’_{i}a_{i}=1$,当i>1时,$\bf a’_{i}a_{j}=0(i=1,2,\dots,i-1)$
(2)$Var(Z_{i})=\mathop{Max{Var(a’\bf X)}}\limits_{\bf a’a=1,\bf a’a_{j}=0(j=1,2,\dots,i-1)}^{}$
刚称$Z_{i}$是原始指标的第$i$主成分($i=2,\dots,m$)。
由上述定义可知,当$i\neq j$时,主成分$Z_{i}$与$Z_{j}$是互不相关的,并且$Z_{1}$是原始指标$X_{1},X_{2},\dots,X_{m}$的一切线性组合中方差最大者,$Z_{2}$是与$Z_{1}$不相关的、除$Z_{1}$以外的$X_{1},X_{2},\dots,X_{m}$一切线性组合中方差最大者,$Z_{m}$是与$Z_{1},Z_{2},\dots,Z_{m}-1$都不相关的,除$Z_{1},Z_{2},\dots,Z_{m}-1$以外的$X_{1},X_{2},\dots,X_{m}$一切线性组合中方差最大者。从理论上讲,求得的主成分个数最多可能有m个,这时,m个主成分就反映了全部原始指标所提供的信息,鉴于主成分分析的目的主要是用较少的综合指标来反映全部原始指标中的主要信息,因此在实际工作中,所确定的主成分个数总是小于原始指标的个数。
为方便讨论,以m=2为例说明主成分分析的几何意义,设个体具有两个观测指标X1和X2,它们之间具有较强的相关性,测量n例这样的个体的值,将所得的n对数据在以X1为横轴,X2为纵轴的二维坐标平面中的苫,得到如下的散点图:
由上图可以看出,由于$X_{1}$和$X_{2}$具有较强的相关性,这n个点的分页呈现出直接化的趋势;同时它们沿$X_{1}$轴方向和$X_{2}$轴方向都具有较大的变异变。个体在某个方向上的变异度可以用该方向上相应观测变量的方差来定量地表示。显然,如果只考虑$X_{1}$和$X_{2}$中任何一个方向上的方差,就将损失原始观测数据中很大一部分信息。如果我们将坐标轴$X_{1}$和$X_{2}$同时按逆时针方向作一个放置,得到新的坐标轴$Z_{1}$和$Z_{2}$,使得在亲折坐标平面上,这n个点的分布基本上不再具有有相关性,且它们的变异主要集中在Z1方向上,而在Z2方向上的变异较小,此时若取$Z_{1}$作为第一主成分,则$Z_{1}$就反应了原始指标$X_{1}$和$X_{2}$所包含的主要信息。
由主成分的定义可知,各主成分互不相关,即任意两个主成分$Z_{i},Z_{j}$的协方差为0,即
且各主成分的方差满足:
于是由公式(1-2)定义的随机向量$Z$的协方差矩阵为:
由主成分定义中的条件(1)可知,这里的方阵$A$是正交阵,即$A’A=I$(I为单位矩阵),由此可得:
由上述公式可知,求原始指标$X_{1},X_{2},\dots,X_{m}$的主成分问题,实际上就是要求满足上述条件的正交阵$A$,即随机微量$X=(X_{1},X_{2},\dots,X_{m}$的协方差矩阵$Cov(X)$的特征值(eigenvalue)与特征向量(eigenvector)。
下面讨论怎么由一组$X_{1},X_{2},\dots,X_{m}$的样本观测值求出主成分,假设收集到的原始数据共有n例,每例测得m个指标的数值,记录如下所示:
将原始指标标准化(标准化通俗讲就是每一列的每个数字每去这一列的均值,然后除以标准差),然后用标准化的数据$X’_{ij}$来计算主成分,为了方便计算,下面的公式中仍用$X_{ij}$表示标准化后的指标数据,$\bf X$为标准化后的数据矩阵,则:
标准化后,$X$的相关矩阵即为协方差矩阵$Cov({\bf X})$
对角线上分别是$r_{11}$,$r_{22}$和$r_{mm}$的方差,非对角线上是协方差。协方差是衡量两个变量同时变化的变化程度。以两个变量x、y来举例,例如协方差大于0,表示x和y若一个增,另一个也增;小于0表示一个增,一个减。如果x和y是统计独立的,那么二者之间的协方差就是0;但是协方差是0,并不能说明x和y是独立的。协方差绝对值越大,两者对彼此的影响越大,反之越小。协方差是没有单位的量,因此,如果同样的两个变量所采用的量纲发生变化,它们的协方差也会产生数字上的变化。
由公式1-3得知,求主成分的问题,实际上是求出$X$的协方差矩阵$Cov(X)$(这里即为$X$的相关矩阵$R$)的特征值和特征向量,由于$R$为半正定矩阵,故可由R的特征方程
解得每一特征值$\lambda_{i}$对应的单位特征向量$a_{i}=(a_{i1} a_{i2} \dots a_{im}’$,从而求得各主成分,即
有关特征值和特征向量的理解可以参考后文引用的文章。这里简单提一下,如果把矩阵看作是运动,对于运动而言,最重要的就是运动的速度和方向,那么特征值就是运动的速度,特征向量就是运动的方向。
即$Z_{i}$与$Z_{j}$的相关系数为0,即:
因此各主成分间的相关系统矩阵为单位矩阵。
可以证明,各原始指标$X_{1},X_{2},\dots,X_{m}$的方差和与各主成分$Z_{1},Z_{2},\dots,Z_{m}$的方差和相等,将数据标准化后,原始指标的方差和为$\sum\limits_{i=1}^m \lambda_{i}$,即有$m=\sum\limits_{i=1}^m \lambda_{i}$
各指标所提供的信息量是用其方差来衡量的。由此可知,主成分分析是把m个原始指标$X_{1},X_{2},\dots,X_{m}$的总方差分解为m个互不相关的综合指标$Z_{1},Z_{2},\dots,Z_{m}$的方差之和,使第一主成分的方差达到最大,即变化最大的方向微量所相应的线性函数,最大方差为$\lambda$。其中${\lambda}/{\sum\limits_{i=1}^m \lambda_{i}}$表明了第一主成分$Z_{1}$的方差在全部方差中所占的比值,称为第一主成分的贡献率,这个值越大,表明$Z_{1}$这个指标综合原指标$X_{1},X_{2},\dots,X_{m}$的能力越强,也可以说,由$Z_{1}$的差异来解释$X_{1},X_{2},\dots,X_{m}$的差异的能力越强,正是因为这一点,才把$Z_{1}$称为$X_{1},X_{2},\dots,X_{m}$的第一主成分,也就是$X_{1},X_{2},\dots,X_{m}$的主要部分,了解到这一点,就可以知道为什么主成分是按特征值$\lambda_{1},\lambda_{2},\dots,\lambda_{m}$进行排序的。
通常情况下
为第i主成分的贡献率;而称
为前k个主成分的累积贡献率。
通常并不需要全部的主成分,只用其中的前几个,一般来说,主成分的保留个数按以下的原则来进行:
当前k个主成分的累积贡献率达到某一特定值时(一般以大于70%为宜,有的时候会要求大于80%),则保留前k个主成分。
即若主成分$Z_{i}$的特征值$\lambda_{i} \geq 1$,则保留$Z_{i}$,否则就去掉该主成分。这个与碎石图类似。
为了了解各主成分与各原始指标之间的关系,在主成的表达式中,第$i$主成分$Z_{i}$的特征值的平方根$\sqrt{\lambda_{i}}$与第$j$原始指标$X_{j}$的系数$a_{ij}$的乘积,即
为因子载荷(factor loading),由因子载荷构成的矩阵称为因子载荷阵,事实上因子载荷$q_{ij}$就是第$i$主成分$Z_{i}$与第$j$原始指标$X_{j}$之间的相关系数,它反映了主成分$Z_{i}$与原始指标X_{i}$之间联系的密切程度与作用的方向。
对于具有原始指标测定值$(X_{i1},X_{i2},\dots,X_{im})$的任一样品,可先用标准化变换式将原始数据标化,即:
然后代入各主成分的表达式,即
求出该样本各主成分值,这样求得的主成分值称为该样本的主成分得分,利用样品的主成分得分,可以对样品的特征进行推断和评价。
例22-1 某研究者测得84名10岁男孩的身高、坐高、体重、胸围、肩宽、肺活量等6项生长发育指标,数据见表22-2。试作主成分分析。
|
|
剔除掉第1列后的数据:
|
|
|
|
在上图中,其中第二行,即Proportion of Variance表示的是贡献率;而第3行,即Cumulative Proportion表示的是累积贡献率;由这个表可以看出,主成分取3个比较合适,此时的累积贡献率为88.92%,接近于90.00%。
|
|
这张表的结果可以得出各个主成分的与相应变量的系数,如果取3个主成分,则可以得出前三个主成分为:
前面的第1到第5步用的是R中的base的函数,没有找到计算因载荷的矩阵方法,现在用psych包来进行计算。
|
|
psych包中的fa.parallel用来提取主成分,绘制的图形如下所示:
这种图叫碎石图,把对应各个主成分的特征值按从大到小的顺序在图上绘出,选取主成分个数至碎石图发生斜率明显变化为止。其中,蓝线是基于观测特征值的碎石检验,根据100个随机数据矩阵推导出来的特征值均值是红线,从图中可以看出:碎石检验图形最大变化处上面只有一个成分;特征值大于随机模拟数据的也只有一个主成分;特征值大于1的也只有一个主成分。所以这一组数据使用一个主成分即可保留数据集的大部分信息(这一点与书上的不符,但主成分的挑选我觉得参入了一定的主观成分在里面,选取几个主成分,要综合几个因素来看,书上选取的是3个主成分,为了保持一致,我们后面的分析也选3个主成分)。
提取3个主成分,如下所示:
结果如下所示:
红色框中就是因子载荷矩阵。
由因子载荷矩阵可知,第一主成分$Z_{1}$在各原始指标上的因子载荷较为均匀,故可认为该主成分反映的是各原始指标的综合信息;第二主成分$Z_{2}$在$X_{1}$(身高)、$X_{3}$(坐高)及$X_{4}$(胸围)上的因子载荷较大,故可认为该主成分反映的是体型方面的信息;而第三主成分$Z_{3}$则主要反映了来自原始指标$X_{6}$(肺活量)的信息。此外,还可知道,第一主成分$Z_{1}$与各原始变量之间的关系较为密切,第二主成分$Z_{2}$与原始变量$X_{1}$、$X_{2}$及$X_{4}$之间的关系较为密切,而第三主成分$Z_{3}$与原始变量$X_{6}$之间的关系较为密切。
此外,结果中还有h2和u2以及com,它们的解释如下:
h2栏指成分公因子方差,主成分对每个变量的方差解释度;
u2栏指成分唯一性,就是方差不能被主成分解释的比例。显然u2=1-h2;
当我们提取的主成分不止一个时,使用主成分旋转会使各主成分所代表的实际意义更容易被解释。主要的旋转方法包括以下两种:
|
|
经过主成分旋转的三个主成分仍然不相关,对变量的解释性不变,这是因为变量的群组没有发 生变化。另外,三个主成分旋转后的累积方差解释性没有变化(89%),变的只是各个主成分对方差的解释度(第一主成分从71%变为39%,第二主成分从10%变为30%,第三主成分从8%变为20%)。
获取主成分与原始变量的线性关系(保存在模型的weights部分):
|
|
这样,我们就可以得到主成分与原始变量之间的相关关系:
RC1 = -0.18X1 - 0.273X2 + 0.47X3 + …
RC2 =0.535X1 + 0859X2 - 0.067X3 + …
RC3 = 0.049X1 - 0.267X2 - 0.173X3 + …
上面的代码只是为了说明PCA的原理,在实际绘图中可以使用其他的R包来进行PCA图的生成,下面是一段生成PCA的代码:
|
|
这是对沈梦圆PPT的整理。是一篇非常好的关于RNA-seq数据分析的文章。
mRNA在生物个体内RNA的组分中只占很小的一部分,rRNA占绝大多数。一般我们说RNA-seq指的都是mRNA-seq,后面的流程也都是主要针对mRNA-seq数据分析的。在科学家们的努力下,可以把那些非编码RNA提取出来建库,进行测序。
一个成功的RNA-seq研究,起决定性因素的是一个好的实验设计。还依赖于建库的类型、测序深度和设置适于的生物重复。并且尽量减少测序本身以外带来的数据误差。
如下所示:
需要注意的问题:
1.一般生物体中的的RNA中,rRNA占绝大多数,含量超过90%,而mRNA的含量在1-2%左右。对于真核生物,一般使用加poly(A)选择性富集mRNA或者而原核生物则是通过去除rRNA;
2.是否建stand-preserving库;
3.对于Illumina,测序插入片段一般小于500bp。确定合适长度的插入片段是后续测序和分析的关键;
4.单端还是双端测序毫无疑问的是,单端测序更便宜一些,如果你研究的某个物种的基因表达水平,并且它的转录组已经被注释很好了,单端测序产生的数据量一般是足够的了。双端测序呢,它的读长更长,更适合于那些没有被注释的转录组物种的研究,便于其转录本的从头拼接。
误差分为技术误差和生物学差异。技术误差-可以通过选择最优化的实验测序程序;生物误差-三个生物学重复。
然后设定生物学重复对差异基因的检出率(真阳性率,TPR)的提高具有明显效果。上面说增加测序深度可以检测到低丰度基因,但是对任何样品来说的当测序深度增加再增加,它就会到达平台期。由于科研经费有限,无法无限制地增加样本数或数据量。
所以在生物学重复数和单个样本测序量上必须找到平衡点。在总数据量不变的情况下,将总数据量分配到更多的生物学重复样本中,差异分析结果的可靠性在不断提升。对于RNA-seq,生物学重复数的价值要大于单个样本测序量。但增加生物学重复的样本数,意味着要增加建库费用。因此,即使总数据不变,设置过多的生物学重复也是不合理的。
我们最终确定设置多少生物学重复还是需要看样本个体之间的差异大不大,这点我们一般都很清楚,在测序之前,如果你所研究的现象在两个实验样本之间差异很稳定的话,就可以少设置一些重复,差异不稳定的话有时候设置10个/20个都不够。具体问题具体分析.
测序深度(Sequencing depth),也叫乘数,指每个碱基被测序的平均次数,是用来衡量测序量的首要参数。研究表明,增加测序深度,测序量从1.6M条reads增加到20M条reads,(75bp)但到10M条reads时就已经达到平衡了,80%的基因转录本被检测到。在此基础上增加测序量,它们会比对到已经存在的转录本上。因此即使提高测序深度,低表达水平的基因的检测是比较困难的。并且提高测序深度确实能够增加基因差异表达的敏感度,但是并不能保证检测到的差异具有生物学意义。
直观一些说,如果某个基因在RNA-seq结果显示差异表达,但QPCR结果表明这个基因表达差异不显著,可以认为这个基因RNA-seq结果为假阳性;反之,这个结果就是真阳性。
生物学重复对差异表达分析的影响
在单样本测序量保持不变的情况下,随着生物学重复(n)的提高,差异分析的假阳性率(FP R)基本稳定,但真阳性率(TPR)在不断提高。也就是说提高生物学重复数,实验对差异基因的检测更加敏感,那些差异倍数较小或差异量较低的差异表达基因(此类基因的差异检测难度较大)能够更加容易被检测到。
判定差异分析结果可靠性的指标
假阳性与真阳性
假阳性率(FPR):真实非差异表达的基因中,被错误判定为差异表达的比例,FPR越低越好。
真阳性率(TRP):真实差异表达的基因中,能够正确判定为差异表达的比例,TRP越高越好。
提高生物学重复数,实验对差异基因的检验更加敏感。那些差异倍数较小或差异量较低的差异表达基因(此类基因的差异检测难度较大)能够更加容易被检测到。
在一定的生物学重复数(n)的情况下,随着单样本测序量(Depth)的提高(25% → 100%),真阳性率(TPR)都只有有限的提高(图1)。例如在n=3的情况下,单个样本的测序量从25%提高到100%,TPR仅仅从6.24%提高到8.95%。在表3中,如果Depth等于25%不变,当n从2提高到12,TPR的提高则是非常明显的。因此测序深度对结果改善效果不如增加生物学重复。
总数据量不变,生物学重复数与单样品测序量最佳组合
如果保持总测序量不变(即如果生物量重复数为n,则单个样品的测序量降低为1/n,总数据量为n*1/n=1,保持不变)。如图A,灰色实线代表不同的生物学重复数(n)和单样本数据量(1/n)组合的情况下,真阳性率(TPR)的变化。结果表明,随着n的提高,TPR率不断提高。例如n=2,TPR约为3%,如果n=6,TPR则提高到24.3%。
同时我们也可以对“单样本测序量对差异表达分析的影响”再进行深入观察。如果n保持不变,但单个样本的数据量不断降低,TPR的降低十分缓慢。例如,n=3,单个样本的数据量从100%降低到15%,TPR的值一直处于平台期,仅仅从9%降低到5%。 但是不同的生物学重复数和单样本测序量的组合,对假阳性率( FPR)的影响却较小。如图B,灰色实线代表不同生物学重复数(n)和单样本数据量(1/n)组合的情况下,真阳性率(FPR)的变化。虽然 n 从2 变化到 96,FPR 基本没有太大变化。
从图中我们很容易发现,基于负二项分布的差异分析检验(P value),FPR 对生物学重复数和单个样本数据量均不敏感,始终保持低于 0.1%水平。或者说,这个算法对 FPR 的控制还是非常理想的。
随着测序单价的下降,目前市场上 RNA-seq 类项目的单样本测序量正在不断提高。以 2G,PE100 测序的表达谱项目为例,其对应的测序量为 20M 条 reads。如果一条长度为 1kbp 的低表达基因的表达量为 RPKM=0.5,其理论上可以检测到的 reads 数为 20×0.5=10。所以低丰度基因的检测,对 RNA-seq 这个技术来说并非最大问题。
第二个问题“转录本表达量的高低变化”比“转录本的有无”更具有普遍的生物学意义。虽然个别基因的表达量变化程度,可以使用 Qpcr 来验证。但我们往往也使用所有差异基因来统计某些规律。例如使用差异基因的 pathway 富集分析来寻找与性状相关的 pathway。如果在全局水平的差异基因集并不可靠,那么 pathway富集分析得出的结论的可靠性自然也受到影响。而全局水平的差异基因数量巨大,是难以使用 Qpcr 验证的。因此,定量以及差异分析的准确性是在 RNA-seq 中更值得关心的问题。
RNA-seq文库的制备和测序过程:RNA碎裂,cDNA合成,接头连接,PCR扩增,加标签(多样品混合测序),上泳池测序;
原始数据回来后,你做完备份以后,做的第一件事情就是看看数据质量如何,一般来自llumina测序平台用软件FastQC看;其他平台的数据用软件NGSQC。一般会有原始数据的序列质量,GC含量,存在的接头以及K-mers子串图并且重复序列太多的reads。并且reads 3‘末端的质量低于前段,原因是随着测序读长的增加,酶活性下降,荧光强度也在下降,因此测序数据质量逐渐降低乃是自然趋势。常用的数据过滤的软件有FASTX-Toolkit and Trimmomatic,其他还有许多,你也可以自己写代码处理数据。
比对上的reads占总reads的百分比; Reads比对到外显子和参考链上的覆盖度是否一致;比对到基因组序列:多重比对reads?比对到转录组序列:来自未被注释的转录本的reads会丢失; 产生更多的多重比对reads; 转录本被定量以后,应该看一下GC含量和基因长度偏差,确定定量的方法是否适用。
需要注意的问题:
把所有样本的reads混合用于转录本的拼接。二代测序的转录组reads用于拼接还是存在一些问题的,最终拼接结果不太理想。一个转录本的拼接结果会是10~100contigs。三代测序的读长直接可以把一个转录本读完了,完全不需要拼接。
RPKM/FPKM/TPM用来表示RNA-seq基因表达水平的值;对于单端测序RPKM和FPKM值是一样的,FPKM可以转换成TPM。Cufflinks(支持双端测序数据,并且需要GTF格式的注释文件)。
功能分析是标准转录组分析流程的最后一步,分析差异表达基因的分子功能和代谢通路。
参考资料: