一个小任务,证明在某个文献看到的这句话: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的。意思是希望大家明白,可观规律是很难因为数据库版本更新而改变的! 其中浙江大学的学徒贾方完成度最好,所以分享一下他八次尝试完成作业的学习历程!
cat gencode.v32.chr_patch_hapl_scaff.annotation.gtf|awk '{if($3=="transcript")print $5-$4}'>trc_length.txt
transcript
的行取出来,然后坐标相减。出来的结果没啥问题,最开始我就是使用这样简单粗暴的坐标之差来定义转录本长度!
$ 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 #最大值
首先怀疑是我的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 #中位数
$ 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 #最大值
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).
$ 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 #中位数
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
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))
transcriptLengths
这个函数计算得到的却是不一样的。transcriptLengths
这个函数的算法不是单纯的非冗余外显子之和,但是帮助文档里说:The length of a processed transcript is just the sum of the lengths of its exons.
这里很令人疑惑。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
> median(trx_len)
[1] 868
> max(trx_len)
[1] 205011
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
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
length(unique(tmp$trx)) [1] 145126 median(trx_len) [1] 909 max(trx_len) [1] 108861 ```
length(unique(tmp$trx)) [1] 8123 median(trx_len) [1] 2388 max(trx_len) [1] 102741 ```