第八次尝试，他终于发现了人类参考基因组转录本的秘密
最近在学徒群给了他们一个小任务，证明我在某个文献看到的这句话：The median length of human transcripts is 2186 nt, with the longest transcripts having sizes of up to 101,206 nt. (These numbers
are based on UCSC hg19 annotation.)
我希望他们可以基于gencode的v32也测试看看，如果不行，再去找hg19的。意思是希望他们明白，可观规律是很难因为数据库版本更新而改变的！其中浙江大学的学徒贾方完成度最好，所以分享一下他八次尝试完成作业的学习历程！
第一次尝试：基于gencode v32的gtf文件
step1 下载gtf，提取长度信息
在gencode官网下载了v32的gtf文件，解压后用下面的代码得到所有transcript的长度信息。
cat gencode.v32.chr_patch_hapl_scaff.annotation.gtf|awk '{if($3=="transcript")print $5-$4}'>trc_length.txt
计算方法是：把注释为
transcript的行取出来，然后坐标相减。检验了一下前几个，没有问题（前20列共3个transcript，且14409-11869=2540，13670-12010=1660，29570-14404=15166）
出来的结果没啥问题，最开始我就是使用这样简单粗暴的坐标之差来定义转录本长度！
step2 找到中位数和最大值
$ wc trc_length.txt
248592 248592 1357637 trc_length.txt
$ cat trc_length.txt |sort -n|awk 'NR==124296{print $0}'
10259 #中位数
$ cat trc_length.txt |sort -n|awk 'NR==248592{print $0}'
2471656 #最大值
总共24万多个转录本，中位数和最大值都差得非常远。
step3 试图纠错
首先怀疑是我的gtf文件有问题
刚才用的是中间这个，现在用另外两个进行同样的操作
$ wc trc_length_chr.txt
227462 227462 1248010 trc_length_chr.txt
$ cat trc_length_chr.txt |sort -n|awk 'NR==113731{print $0}'
11030 #中位数
$ cat trc_length_chr.txt |sort -n|awk 'NR==227462{print $0}'
2471656
$ wc trc_length_pri.txt 227529 227529 1248322 trc_length_pri.txt
$ cat trc_length_pri.txt |sort -n|awk 'NR==227529{print $0}'
2471656 #最大值
$ cat trc_length_pri.txt |sort -n|awk 'NR==113764{print $0}'
11026 #中位数
$ cat trc_length_pri.txt |sort -n|awk 'NR==113765{print $0}'
11026 #中位数
结果差别不大，看来不是这个问题，继续怀疑是需要区分protein-coding与否。
$ cat gencode.v32.annotation.gtf|tr '"' 'a'|awk '{if($3=="transcript" && $18=="aprotein_codinga;")print $5-$4}'>trc_length.txt #由于protein_coding这个信息在原文件中加了双引号："protein_coding";而awk匹配不到，所以这里用了一个笨办法，干脆把所有双引号替换成字符a，这样就能匹配到了。
$ wc trc_length.txt
83986 83986 489050 trc_length.tx
$ cat trc_length.txt |sort -n|awk 'NR==41993{print $0}'
20997 #中位数
$ cat trc_length.txt |sort -n|awk 'NR==83986{print $0}'
2471656 #最大值
但是答案还是不对。
第二次尝试：基于ensembl的gtf文件
本来想用原文提到的这个UCSC hg19 annotation再试一下，不过搜索这几个关键词的时候看到biostars上一个老兄这么一番话：
Is there a reason you want to use the UCSC annotation? The one from Ensembl/Gencode is almost always better (there's a reason that UCSC now uses the copy from gencode).那还是用ensembl试一下吧，还是一样的代码。
$ cat Homo_sapiens.GRCh38.90.chr.gtf|awk '{if($3=="transcript")print $5-$4}'>trc_length_ensembl.txt
$ wc trc_length_ensembl.txt
200243 200243 1086390 trc_length_ensembl.txt
$ cat trc_length_ensembl.txt |sort -n|awk 'NR==200243{print $0}'
2471656 #最大值
$ cat trc_length_ensembl.txt |sort -n|awk 'NR==100121{print $0}'
9634 #中位数
$ cat trc_length_ensembl.txt |sort -n|awk 'NR==100122{print $0}'
9635 #中位数
果然依旧是相似的结果。
第三次尝试：基于R包：`TxDb.Hsapiens.UCSC.hg19.knownGene`
在汇报了上面的结果之后，曾老师提示我：转录本长度，不是坐标相减去，而是外显子之和。
于是我想到这个和基因长度有点像，之前做RNA-seq流程的时候遇到过基因长度的问题，当时是没有完全解决的，用公司给的信息凑合了，当时还查到生信菜鸟团的这篇教程（http://www.bio-info-trainee.com/3298.html），是用perl命令计算非冗余CDS之和。现在我的需求是从gtf文件计算出转录本长度（对应的外显子长度之和），以我目前的shell命令水平较难解决，因此想用教程中提到的R包解决这个问题。
参考：http://www.bio-info-trainee.com/3991.html；https://mp.weixin.qq.com/s/l7fTh9MMFCXLWMY95eSjvg
最直接的方法是用
transcriptLengths这个函数
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
t_l=transcriptLengths(txdb)
t_l=na.omit(t_l)
tmp=sort(t_l$tx_len)
> median(tmp)
[1] 2377
> max(tmp)
[1] 109223
结果和原文类似，但不完全一致，而且这种方法有一个问题，就是只得到了73432个转录本（去除na值之前也只有82960个），所以这个包里所谓的转录本信息应该是没有取全部转录本，而是只取了protein-coding的转录本。
以下代码是我想用类似计算每个基因的非冗余外显子长度之和的方式计算每个转录本对应的非冗余外显子长度之和。
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
exon_txdb=exons(txdb)
trx_txdb=transcripts(txdb)
o = findOverlaps(exon_txdb,trx_txdb)
t1=exon_txdb[queryHits(o)]
t2=trx_txdb[subjectHits(o)]
t1=as.data.frame(t1)
# mcols作用是提取 a DataFrame object containing the metadata columns
t1$trxid=mcols(t2)[,1]
t = lapply(split(t1,t1$trxid),function(x){
tmp=apply(x,1,function(y){
y[2]:y[3]
})
length(unique(unlist(tmp)))
})
tmp=data.frame(trx_id=names(t),length=as.numeric(t))
最后的结果和之前的方法对比，发现是有问题的，比如编号为1，2，3的这三条转录本，它们来自同一个基因，因此对应同样的3个外显子，但是用我的方法计算得到的长度是一样的，而直接用
transcriptLengths这个函数计算得到的却是不一样的。
这就说明
transcriptLengths这个函数的算法不是单纯的非冗余外显子之和，但是帮助文档里说：The length of a processed transcript is just the sum of the lengths of its exons.这里很令人疑惑。但是最大的问题还是，用这个R包终究只能得到82960个转录本的长度，所以正道应该还是用gtf文件。
第四次尝试：基于gencode v32的gtf文件
将上面的内容汇报后，曾老师给的提示是
我在这些题目里看到gtf是可以载入R操作了，之前我居然没想到，于是载入R尝试实现刚才的需求：从gtf文件计算出转录本长度（对应的外显子长度之和）。
这里用的gtf文件时CHR的这个，因为了解了一下
代码：
gtf <- read.table('gencode.v32.annotation.gtf.gz',stringsAsFactors = F,
header = F,comment.char = "#",sep = ' ')
exon <- gtf[gtf[,3]=='exon',] #首先筛选到了exon信息
exon$trx <- lapply(exon[,9], function(x){
y=strsplit(x,';')[[1]][2]
strsplit(y,' ')[[1]][3]
}) #把gtf的第9列拆一下，留下转录本编号作为exon的分组信息
exon <- exon[,c(4,5,10)]
head(exon)
## V4 V5 trx
## 3 11869 12227 ENST00000456328.2
## 4 12613 12721 ENST00000456328.2
## 5 13221 14409 ENST00000456328.2
## 7 12010 12057 ENST00000450305.2
## 8 12179 12227 ENST00000450305.2
## 9 12613 12697 ENST00000450305.2
length(unique(exon$trx))
## [1] 227462
# 20多万个转录本
exon$trx=unlist(exon$trx)
exon$trx=as.factor(exon$trx) #为后续用tapply做准备
exon$len=exon[,2]-exon[,1]
trx_len=tapply(exon[,4],exon[,3],sum) #分组计算和，后来想到这里也可以用group_by + summarise
median(trx_len)
## [1] 876
max(trx_len)
## [1] 205011
可惜，差得还是很多。
第五次尝试：基于gencode v32 gtf文件的GRCh37兼容版本
汇报后经前辈提醒，注意到应该使用hg19基因组，因此特意重下了gencode v32 gtf文件的GRCh37兼容版本，重新跑了一遍上面的R代码。
结果是
> median(trx_len)
[1] 868
> max(trx_len)
[1] 205011
和GRCh38几乎没有差别。
第六次尝试：基于转录本fastq文件
汇报后曾老师的提示：gencode，ensembl，ccds数据库的转录本fa文件，可以直接统计碱基数量，跨过gtf问题。
因此我下载了gencode的转录本fa文件来尝试，还是用R来完成，代码是（其实读取fa文件应该是使用专门的R包，我这个方法太拙劣了）
rm(list=ls())
options(stringsAsFactors = F)
trx <- read.table('gencode.v32lift37.transcripts.fa',stringsAsFactors = F,
header = F,comment.char = ">",sep = '\n')
len <- apply(trx,1,nchar)
#本来以为这时候就大功告成了，但我发现长度几乎都是60，说明fasta文件同一个转录本的序列被R识别成了多行，每行60个碱基。
#思来想去决定用for循环解决这个问题，找到不是60的元素时，把它作为转录本之间的分界，把这之前的数全部加起来，当然这样做默认转录本长度不会是60的倍数。
# tmp=head(len,83)
# table(tmp==60)
# count=rep(0,5)
# index=1
# for(i in 1:length(tmp)){
# if(tmp[i]==60) count[index]=count[index]+60
# if(tmp[i]<60){
# count[index]=count[index]+tmp[i]
# index=index+1
# }
# }
table(len==60)
count=rep(0,225274)
index=1
for(i in 1:length(len)){
if(len[i]==60) count[index]=count[index]+60
if(len[i]<60){
count[index]=count[index]+len[i]
index=index+1
}
}
#还好，由于循环体足够简单，这个600多万次的循环1s不到就计算出来了。
median(count)
max(count)
可惜最后的答案依旧是
> median(count)
[1] 886
> max(count)
[1] 205012
这个时候由于我用好几种方法计算的结果都是差不多的，我开始有点怀疑原来那句话的真实性。
第七次尝试：取转录本子集
区分protein-coding与否
经曾老师提醒，准备尝试一下区分是否编码蛋白这个特征
rm(list=ls())
gtf <- read.table('gencode.v32lift37.annotation.gtf',stringsAsFactors = F,header = F,comment.char = "#",sep = ' ')
exon <- gtf[gtf[,3]=='exon',]
# exon[1,9]
# tmp=exon[1:50,]
# ttmp=tmp[grepl("protein_coding",exon[1:50,9]),]
exon <- exon[grepl("protein_coding",exon[,9]),] #和之前的代码几乎一样，只是这里加了一个筛选步骤
exon$trx <- lapply(exon[,9], function(x){
y=strsplit(x,';')[[1]][2]
strsplit(y,' ')[[1]][3]
})
tmp <- exon[,c(4,5,10)]
> length(unique(tmp$trx))
[1] 153249
#剩下15万个转录本
tmp$trx=unlist(tmp$trx)
tmp$trx=as.factor(tmp$trx)
tmp$len=tmp[,2]-tmp[,1]
trx_len=tapply(tmp[,4],tmp[,3],sum)
median(trx_len)
max(trx_len)
> median(trx_len)
[1] 961
> max(trx_len)
[1] 108861
中位数还是和原来差不多，最大值接近了。
第8次尝试：区分HAVANA和ENSEMBL
这个时候我想到之前有个事情一直被我忽略了，就是gtf文件里的信息有两个来源，HAVANA和ENSEMBL，可能需要区分，尝试了一下分别算HAVANA的所有转录本/编码蛋白的转录本以及ENSEMBL的所有转录本/编码蛋白的转录本，代码和上面一样，加了一个筛选，结果如下：
length(unique(tmp$trx))
[1] 145126
median(trx_len)
[1] 909
max(trx_len)
[1] 108861
```
ENSEMBL，all
length(unique(tmp$trx))
[1] 8123
median(trx_len)
[1] 2388
max(trx_len)
[1] 102741
```
这么看下来最后一种最接近。
小结
开始做这个作业的时候我万万没想到一句这么简单的话的证明居然会这么复杂，最后虽然没有得到完全一致的答案，但是也对各类转录本的长度特征有了一些了解，重要的是感觉自己解决问题的能力得到了提升，非常感谢曾老师和另一位前辈在这个过程中给我的帮助，这可能才是真正的授人以渔。
备注：另一位前辈就是九月学徒，杨胖子，有诸多笔记，包括九月学徒转录组学习成果展（3万字总结）（上篇）和这个胖子说要更新100篇单细胞数据处理教程
后记
"We cannot teach people anything; we can only help them discover it within themselves."
~Galileo Galilei
可以选择参加11月23号周六下午3-6点，上海张江高科的义诊，https://mp.weixin.qq.com/s/xIOzD4XT8_G8LCZP0kZxgw
或者11月24号周日上午9-12点的拍卖行，两种形式都可以帮助你解决生物信息学问题，https://mp.weixin.qq.com/s/wNi-Ips8WKffAAZAmpo3NA
同时我们还会提供一些招聘见面会，比如 https://mp.weixin.qq.com/s/0i51Fv5pdu4f8W06RGIdGQ