task5b-验证lncRNA只有部分具有polyA尾结构

题目:

下载人和鼠的gtf文件，以及转录本fasta序列文件，自己去探索一下：

gtf文件记录了多少个基因，多少个是蛋白编码基因;
多少个是lncRNA呢？其中各自的具有polyA尾结构的比例是多少呢？

背景知识

真核生物的mRNA都是有polyA尾巴结构
lncRNA只有部分具有polyA尾结构

数据下载

Genocode刚好有人和小鼠的gtf文件以及转录本fasta序列文件

#human
wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.transcripts.fa.gz'
wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz'
nohup wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_38/gencode.v38.metadata.PolyA_feature.gz' &
#mouse
wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/gencode.vM27.chr_patch_hapl_scaff.annotation.gtf.gz'
wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/gencode.vM27.transcripts.fa.gz'
wget -c 'http://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M27/gencode.vM27.metadata.PolyA_feature.gz'

Genocode中gtf文件注释来源

两个来源:ENSEMBL和HAVANA
havana 的全称叫做 human and vertebrate analysis and annotation, 是sager 研究所的1个团队，致力于为人和其他脊椎动物提供高质量的基因和转录本注释，havana 团队使用自己研发的软件，手工对人，小鼠，斑马鱼，还有其他一些脊椎动物进行基因和转录本的注释，特别是对于转录本的isoform 和假基因（这两种是目前常见的基因组注释软件所欠缺的地方）进行手工注释；
对于转录本注释的结果，要求必须有对应的实验证据支持，包括cDNA, EST, 蛋白序列等，所以havana 注视的转录本可信度非常高；
human

$ zcat gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/' |cut -f 2 | sort |uniq -c236459 ENSEMBL #7%
3172759 HAVANA #93%
#在gencode中下载的gtf 文件中，93%来源于HAVANA，7%来源于ENSEMBL

mouse

$ zcat gencode.vM27.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/' |cut -f 2 | sort |uniq -c140143 ENSEMBL #7.5%
1729538 HAVANA #92.5%

人类gtf探索

#下载的不同的gtf,最后结果会有差异
$ zcat gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene"' | wc -l
67049 #一共67049个gene$ zcat gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene" && /protein_coding/' | wc -l
22757 #protein_coding gene个数$ zcat gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene" && /lncRNA/' | wc -l
17915 #lncRNA个数

human-polyA探索

#带ployA尾巴的转录本id
zcat gencode.v38.metadata.PolyA_feature.gz |cut -f 1 |sort -u >hg38_polyA_trans_id
$ head hg38_polyA_trans_id
ENST00000484859.1
ENST00000484859.1
ENST00000412666.1
ENST00000412666.1
ENST00000452176.2
ENST00000452176.2
ENST00000450696.1
ENST00000450696.1
ENST00000429505.6
ENST00000429505.6

#人类hg38-gtf中lncRNA的gene_id 转录本 --对应关系
zcat gencode.v38.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="transcript" && /lncRNA/' | perl -nle '{/gene_id\s+\"(.*?)\"\;\s+transcript_id\s+\"(.*?)\"\;\s+gene_type\s+\"(.*?)\"\;/;print "$1\t$2\t$3"}' > gene_id_trans_lncRNA
$ head gene_id_trans_lncRNA
ENSG00000243485.5   ENST00000473358.1   lncRNA
ENSG00000243485.5   ENST00000469289.1   lncRNA
ENSG00000237613.2   ENST00000417324.1   lncRNA
ENSG00000237613.2   ENST00000461467.1   lncRNA
ENSG00000238009.6   ENST00000466430.5   lncRNA
ENSG00000238009.6   ENST00000477740.5   lncRNA
ENSG00000238009.6   ENST00000471248.1   lncRNA
ENSG00000238009.6   ENST00000610542.1   lncRNA
ENSG00000238009.6   ENST00000453576.2   lncRNA
ENSG00000239945.1   ENST00000495576.1   lncRNA

#用R语言对这两个文件做一个merge
gene_id_trans_lncRNA = read.table('gene_id_trans_lncRNA',header = F)[,1:2]
colnames(gene_id_trans_lncRNA) = c('gene_id', 'trans_id')hg38_polyA_trans_id = read.table('hg38_polyA_trans_id',header = F)
colnames(hg38_polyA_trans_id) = 'trans_id'
hg38_polyA_trans_id$others = hg38_polyA_trans_id$trans_idnewdata = dplyr::left_join(gene_id_trans_lncRNA, hg38_polyA_trans_id ,by='trans_id')
newdata = na.omit(newdata)lncRNA_sum = length(unique(gene_id_trans_lncRNA$gene_id))
lncRNA_polyA = length(unique(newdata$gene_id))cat(paste0(round((lncRNA_polyA/lncRNA_sum)*100,2),"%"))

计算出来人类的lncRNA中有30.96%带有polyA尾巴

小鼠gtf探索

$ zcat gencode.vM27.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene"' | wc -l
55416$ zcat gencode.vM27.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene" && /protein_coding/' | wc -l
21885$ zcat gencode.vM27.chr_patch_hapl_scaff.annotation.gtf.gz |awk '!/^#/ && $3=="gene" && /lncRNA/' | wc -l
9950

mouse-gtf探索

用计算人类lncRNA中polyA的比例中类似的方法得到小鼠的lncRNA中带polyA尾巴的比例为:15.45%

可视化-人类小鼠gene，蛋白编码基因，lncRNA，以及lncRNA带polyA尾巴的数量对比

library(ggplot2)
abcd <- read.table(file = "task5b_data.txt",header = T)
ggplot(abcd,aes(x=category,y=number,fill=group))+geom_bar(stat="identity",width = .45,position = "dodge")+ylab("total")+labs(title ="Statistics of Human_vs_mouse")+theme(plot.title = element_text(hjust = 0.5))+geom_text(aes(label=number, y=number+0.05), position=position_dodge(0.6), vjust=0)

参考：

都说lncRNA只有部分具有polyA尾结构，请证明
HAVANA 团队简介