发布时间 : 星期一 文章有参考基因组的转录组生物信息分析模板更新完毕开始阅读d145286eee06eff9aef807d9
2 1000003 1000004 TSS TTS 0 CUFF.100 CUFF.100 1 1 3456744 3457165 3457165 3474806 3478178 3474806 + + 00 00 00 225 225 225 2.00000000ENSGALT000000105.00000000ENSGALT00000010(1) event_id: AS事件编号 (2) event_type: AS事件类型 (TSS, TTS, SKIP_{ON,OFF}, XSKIP_{ON,OFF}, MSKIP_{ON,OFF}, XMSKIP_{ON,OFF}, IR_{ON ,OFF}, XIR_{ON,OFF}, AE, XAE) (3) gene_id: cufflink组装结果中的基因编号 (4) chrom: 染色体编号
(5) event_start: AS事件起始位置 (6) event_end: AS事件结束位置
(7) event_signature: AS事件特征 (for TSS, TTS - inside boundary of alternative marginal exon; for *SKIP_ON,the coordinates of the skipped exon(s); for *SKIP_OFF, the coordinates of the enclosing introns; for *IR_ON, the end coordinates of the long,
intron-containing exon; for *IR_OFF, the listing of coordinates of all the exons along the path containing the retained intron; for *AE, the coordinates of the exon variant) (8) strand: 基因正负链信息
(9) fpkm: 此AS类型该基因表达量
(10) ref_id: 此基因在参考注释文件中的编号
5 新转录本预测
将所有测序reads数据的基因组定位结果放到一起，用 Cufflinks 进行组装，然后用Cuffcompare和已知的基因模型进行比较，可以:(1)发现新的未知基因（相对于原有基因注释文件）；(2)发现已知基因新的外显子区域；(3)对已知基因的起始和终止位置进行优化。新基因和新外显子区域预测结果为GTF格式的注释文件。GTF格式的详细说明可参考(http://mblab.wustl.edu/GTF22.html)
表5.1 新转录本结构注释结果
seqname source feature start end score strand frame 1 1 1 1 novelGene exon 18531 19499 . novelGene exon 20813 21813 . novelGene exon 23917 24402 . novelGene exon 25189 26100 . + + + + . . . . attributes gene_id \\gene_id \\gene_id \\gene_id \\(1) seqname：染色体编号
(2) source：来源标签，这里的novelGene指新基因 (3) feature：区域类型，目前我们预测外显子区域 (4) start：起始坐标 (5) end：终止坐标 (6) score：不必关注 (7) strand：正负链信息 (8) frame：不必关注
(9) attributes：属性，包括基因编号、转录本编号等信息
表5.2 已知基因结构优化
Gene_id ENSGALG00000000003 ENSGALG00000000004 ENSGALG00000000011 ENSGALG00000000013 Chromosome Strand 1 Z 6 22 + - - + Original_span 19947199~19967662 34293201~34299523 30928671~30932616 2783575~2787337 Assembled_span 19947199~19967662 34293201~34299554 30928671~30932616 2783575~2787453 (1) Gene_id：原注释文件中基因命名编号 (2) Chromosome：染色体编号 (3) Strand：正负链信息
(4) Original_span：原注释文件中基因起始位置~终止位置
(5) Assembled_span：转录组拼接结果中基因起始位置~终止位置
6 SNP和Indel分析
SNP全称Single Nucleotide Polymorphisms，是指在基因组上由单个核苷酸变异形成的遗传标记，其数量很多，多态性丰富。从理论上来看每一个SNP 位点都可以有4 种不同的变异形式，但实际上发生的只有两种，即转换和颠换，二者之比为1:2。SNP在CG序列上出现最为频繁，而且多是C转换为T，原因是CG中的C常为甲基化的，自发地脱氨后即成为胸腺嘧啶。一般而言，SNP是指变异频率大于1%的单核苷酸变异。Indel(insertion-deletion)是指相对于参考基因组，样本中发生的小片段的插入缺失，该插入缺失可能含一个或多个碱基。 我们通过samtools和picard-tools等工具对比对结果进行染色体坐标排序、去掉重复的reads等处理，最后通过变异检测软件GATK(McKenna et al., 2010)分别进行SNP Calling和Indel Calling，并对原始结果进行过滤，得到如下表形式的分析结果。其中Indel分析结果每列的含义和SNP结果是一致的。
表6 SNP分析结果
#CHROM 1 POS 502 REF A ALT G HS1 .. HS2 .. HT1 11 HT2 .. HW1 00 HW2 .. 1 1 1 563 1213 1316 C A G A G A .. .. 00 .. .. 00 11 11 01 11 .. .. 00 .. 00 .. .. 00 #CHROM：SNP位点所在染色体 POS：SNP位点坐标
REF：参考序列在该位点的基因型 ALT：该位点的其它基因型 other coloums：每个个体该位点的基因型(0 与REF一致；1 与ALT一致；. 缺少数据支持)
7 基因表达水平分析
一个基因表达水平的直接体现就是其转录本的丰度情况，转录本丰度程度越高，则基因表达水平越高。在RNA-seq分析中，我们可以通过定位到基因组区域或基因外显子区的测序序列(reads)的计数来估计基因的表达水平。Reads计数除了与基因的真实表达水平成正比外，还与基因的长度和测序深度成正相关。为了使不同基因、不同实验间估计的基因表达水平具有可比性，人们引入了RPKM的概念，RPKM(Reads Per Kilo bases per Million reads)是每百万reads中来自某一基因每千碱基长度的reads数目。RPKM同时考虑了测序深度和基因长度对reads计数的影响，是目前最为常用的基因表达水平估算方法 (Mortazavi et al., 2008)。
结果文件分别统计了不同表达水平下基因的数量以及单个基因的表达水平。一般情况下，RPKM数值0.1或者1作为判断基因是否表达的阈值，不同的文献所采用的阈值不同。
表7.1 不同表达水平区间的基因数量统计表
RPKM Interval 0-1 1-3 3-15 15-60 >60 HS1 HS2 HT1 HT2 HW1 HW2 11678(44.62%) 11157(42.63%) 11644(44.49%) 11552(44.14%) 11663(44.56%) 11652(44.52%) 3416(13.05%) 3829(14.63%) 3497(13.36%) 3622(13.84%) 3359(12.83%) 3503(13.39%) 6586(25.17%) 6741(25.76%) 6719(25.67%) 6731(25.72%) 6441(24.61%) 6522(24.92%) 3436(13.13%) 3421(13.07%) 3278(12.53%) 3277(12.52%) 3612(13.80%) 3442(13.15%) 1055(4.03%) 1023(3.91%) 1033(3.95%) 989(3.78%) 1096(4.19%) 1052(4.02%) 表7.2 基因表达水平统计表
geneID 000003 HS1 01677 HS2 40486 HT1 85256 HT2 70462 HW1 19692 HW2 09366 ENSGALG000000.1774627480.57509924280.18207677010.32875021680.22426855480.2831623894ENSGALG000009.4448453447.21882264978.93638040289.694005382810.4063844079.5631290234000004 000011 000013 41505 60835 18205 4696 8464 7827 376 6865 0897 1103 6002 7573 8955 413 2031 8064 986 3859 ENSGALG0000014.7023552914.23548457111.32025669512.29436572213.35758690510.060029486ENSGALG000006.7488688417.48873293736.66826144957.41338358787.53637022217.50828610128 RNA-seq整体质量评估
8.1 表达水平的饱和曲线检查
定量饱和曲线检查反映了基因表达水平定量对数据量的要求。表达量越高的基因，就越容易被准确定量；反之，表达量低的基因，需要较大的测序数据量才能被准确定量。
表达水平的饱和曲线的具体算法描述如下：分别对10%、20%、30%??90%的总体测序数据单独进行基因定量分析，并把所有数据条件下得到的基因的表达水平作为最终的数值。用每个百分比条件下求出的单个基因的RPKM数值和最终对应基因的表达水平数值进行比较，如果差异小于15%，则认为这个基因在这个条件下定量是准确的。