可变剪切的翻译分析
分析AS isoform的编码框
可变剪切产生提前终止的密码子,可能会导致对应功能的缺失或者发生分化;如果这个转录本被翻译出来,就可能导致同源基因在功能上发生分化;重复基因间AS的差异可能会导致基因非功能发生进化,也可能保留一些重复基因
安装EMBOSS
编译的时候确实一个libnucleus.so.6库
下载地址ftp://ftp.pbone.net/mirror/ftp5.gwdg.de/pub/opensuse/repositories/home:/joscott/CentOS_CentOS-6/x86_64/EMBOSS-libs-6.6.0-6.2.x86_64.rpm
##使用yum安装这个库
yum localinstall EMBOSS-libs-6.6.0-6.2.x86_64.rpm
##提示缺失库时,接着下载就行
地址 http://rpm.pbone.net/index.php3/stat/4/idpl/54929788/dir/centos_6/com/EMBOSS-libs-6.6.0-6.2.x86_64.rpm.html
预测isoform的ORF
PacBio文件中保存的是DNA的正链信息
参考基因组中保存的则是有义链的信息
使用merge.gtf中的注释信息提取对应的cDNA序列这一步很关键
##提取对应的isoform序列,这里根据比对的得到的gtf文件提取对应的cDNA序列信息
python ~/github/zpliuCode/script/genestruct/getcDNAsequence.py ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta merge.gtf cDNA.fa
##预测ORF
/public/home/software/opt/bio/software/EMBOSS/6.5.7/bin/getorf -find 1 -noreverse -sequence cDNA.fa
与参考基因组转录本的ORF进行比较
python ~/github/zpliuCode/script/genestruct/compareToreference.py -refgtf Ghirsutum_gene_model.gtf -PBgtf merge.gtf -reforf TM1_reference.orf -PBorf mappingTogene.orf -o test.txt
AS会与细胞质的的NMD机制结合,降解包含有提前终止的密码子。
因此我们在这里分析了Annotion 转录本与unAnnotion转录本的表达水平是存在一个差异;接着对转录本的ORF进行了预测
转录本ORF的统计
这一段可以参考这篇文献
A Comprehensive Analysis of Alternative Splicing in Paleopolyploid Maize
统计annotion转录本中ORF的数据与unannotion转录本的ORF数据
表达水平
转录本长度
ORF长度
可变剪切的比例
根据基因转录本的表达水平,筛选出基因的productive转录本
根据PacBio full-read 数目进行筛选
进而根据基因的表达水平进行筛选
根据筛选出的productive 转录本,统计其他转录本与productive转录本在CDS序列上是否保守;
CDS保守的定义:
CDS长度要一致
CDS序列的保守性达到95%以上
##evolution4/transcriptSVs/TM1
##由于用的是全排列的,并且这个文件只存在CDS长度一样的转录本
awk '$3!=$4{print $0}' A2_genome.txt |sort -k1,4|uniq|sed 's/A2\^//g' >A2_PacBio_CDS.txt
可以发现很多基因存在不同的转录形式,但最终还是产生同样功能的蛋白质。说明这些剪接机制影响的是转录本的UTR区域,可能存在其他的调控作用;根据功能保守,但是存在多种剪接弥散的转录本的基因,根据主要表达转录本的read数占基因总的full-read数进行统计,将基因分为那种转录比较保守的基因、转录后调控比较复杂的基因。
比较同源基因间这种基因的转录后调控是否存在差异。
在检测到PacBio转录本的基因中,有多少基因只转录出单个转录本,而又有多少基因转录出多个转录本
基因组
单个转录本的基因
多个转录本的基因
基因数
转录本数
TM1
13361
18606
31967
83392
A2
6046
14167
20213
67113
D5
6326
11911
18237
51964
比较productive转录本CDS的长度和非productive转录本的CDS长度分布;
这两种转录本似乎在长度上没有什么差异
##根据PacBio read和RNA-seq的大小,获取productive转录本
python ~/github/zpliuCode/transcriptSV/getProductiveTranscript.py ../CO31_32_result/evolution4/TM1_PacBio.txt productive_isoform.txt
awk '{OFS="\t";print $1,$2,int($3),$4,$5,$6}' productive_isoform.txt >1
mv 1 productive_isoform.txt
##获取非保守的转录本
awk '{OFS="\t";print $1,$2,$6,$5,$3}' ../CO31_32_result/evolution4/TM1_PacBio.txt |cat - productive_isoform.txt|cut -f1-5|sort -k1,2 |awk '{OFS="\t";print $3,$4,$5,$1,$2}' |uniq -f3 -u|awk '{OFS="\t";print $4,$5,$1,$2,$3,"noproductive"}'
两种转录本的数目和剪接位点是否是注释的
基因组
Productive
非productive
TM1
31967
51425
A2
20213
46900
D5
18237
33727
分析两种转录本占annotion转录本和非annotion转录本的比例
在鉴定到的productive 转录本中存在0.58-0.83转录本是Annotion的转录本,
而仍旧存在许多未注释的productive转录本,其中D5的比例最少、A2最多;
这可能和基因组的组装质量有关。
基因组
productive
非productive
TM1
22342|9625(0.69)
23785|27640
A2
11780|8433(0.58)
16764|30136
D5
15272|2965(0.83)
19557|14170
行使功能的转录本与其他转录本在CDS长度上的差异
productive转录本CDS保守的转录本也被认为是有功能的转录本,
和功能不一致的转录本的在CDS长度上的差异;
功能一致的转录本与功能发生分化的转录本相比CDS长度明显更长一些;
并且功能一致的转录本也存在差异的剪接方式,说明这些差异的剪接影响UTR区域,对应的是基因的转录调控。分析转录调控在多倍化过程中的差异和同一个基因组内的比例。
大量的功能与主要转录本一致的转录本在 注释信息上是否有差异。(cDNA长度上的差异)气泡图(差异的碱基数映射为点的大小)
##获取功能一致的非productive 转录本
python 01productiveIsoform/01function_same/AddAnnotionTag.py productive_isoform.txt TM1_PacBio_CDS_conserved.txt function_same.txt
##提取这些转录本的完整信息
python ../../CO31_32_result/ORF/AddAnnotionTag.py ../noproductive_isoform.txt function_same.txt 11
cat 11 ../productive_isoform.txt |awk '{print $0"\tfunctionSame"}' >function_same.txt
##与productive 发生功能分化的转录本
cat 11 ../noproductive_isoform.txt |sort |uniq -u|awk '{print $0"\tfunctionDiver"}' >function_Diver.txt
两种转录本的数目和剪接位点是否是注释的
多少比例的基因存在两条以上功能相同的转录本;UTR区域的调控作用
A2 8287个
D5 6780个
TM1 10135个
筛选存在function same转录本前5%的基因进行GO富集分析。
基因组
存在function 的转录本
非function same
TM1
49595(59.47%)
33797
A2
35242(52.51%)
31871
D5
29571(56.90%)
22393
三个基因组中有多少基因存在UTR区域的调控
大约有多少转录本是UTRs 转录本(功能相似但是非productive转录本),比例与拟南芥和水稻类似。
TM1 21.1% (49595-31967)/83392
A2 22.4%
D5 21.8%
大约有27.3%-34.3%的表达的基因存在UTR区域的差异,
分析这些受到调控作用的基因富集在哪些代谢途径中。
基因组
存在UTR差异
不存在
TM1
10135(31.7%)
23238(72.7%)
A2
8287(40.9%)
13278(65.7%)
D5
6780(37.1%)
12576(69%)
剪接位点是否已经被注释
Annotion 转录本数| unAnnotated转录本数
基因组
存在function 的转录本
非function same
TM1
32131|17464
13996|19801
A2
18233|17009
10311|21560
D5
23571|6000
11258|11135
比较productive 转录本与对应与之功能相似转录本的cDNA长度差异。
虽然转录本存在不同的剪接方式,但是仍旧能够产生保守的氨基酸序列,说明基因的剪接影响的是UTR区域(uORF),可能对应着一些转录后的调控,例如miRNA等;影响UTR区域
这个会受测序深度的影响,比如有些非productive 的功能一致的转录本可能会检测不到,导致这个基因不存在UTR区域的调控作用。
获取productive转录本与非productive转录本,但是功能是一致的转录本的cDNA长度
将两个转录本进行比较,看到底是TTS、ployA的差异导致的还是内部的剪接方式差异导致的。
从比较中可以看出两种转录本在剪接方式是否存在差异
##提取两种转录本的gff文件
python 01productiveIsoform/05gffcompared/extract_PacBioAnnotion.py all_isoform.gtf isoformId.txt productive.gtf
##使用gffcompare 进行比较
gffcompare -R -r productive.gtf noproductive.gtf
##统计转录本间的差异
awk '$3~/transc/{print $0}' gffcmp.annotated.gtf |awk -F ";" '{print $(NF-2)}'|sort|uniq -c
splice site剪接位点相同,在TSS和ployA位点上存在差异
splice site 存在差异 ,内部的剪接位点存在差异
UTR发生改变,但是没有影响蛋白质功能;这些基因可能存在UTR区域的调控作用。找几个read数目高的基因或者例子和已经报道的基因。
基因组
splice site完全相同
split site存在差异
TM-1
7970
9627
A2
7725
7241
D5
6255
5066
UTR区域对基因转录调控的调控程度
$1-(a/(a+b))$
a: productive isoform read数目
b:与productive 转录本功能相同但是剪接方式不一致的read
##获取基因中function same转录本的总的full read数目
python 01productiveIsoform/03UTRRegulation/AddAnnotionTag.py
~/work/Alternative/result/Gh_result/CO31_32_result/evolution4/TM1_PacBio.txt ../01function_same/function_same.txt function_same_message.txt
##将基因中功能相似的转录本read相加
awk '{a[$1]+=$6}END{for(i in a){print i"\t"a[i]}}' function_same_message.txt |sort -k1,1 >gene_function_same_read_count.txt
##求每个基因受到UTR区域调控的程度
cut -f1,3 ../productive_isoform.txt|cat - gene_function_same_read_count.txt |sort -k1,1 -k2,2n|awk 'NR%2==0{print $0}NR%2!=0{printf $0"\t"}'|awk '{print $1"\t"$2"\t"$4"\t"1-$2/$4}' >gene_reulate_byUTR.txt
##统计每个调控程度区间的基因数
awk '{for(i=0;i<1;i+=0.1){if($4>=i&&$4<i+0.1){a[i]+=1}}}END{for(i in a){print i"\t"a[i]}}' gene_reulate_byUTR.txt |awk '{print $0"\tA2"}' >plot_gene_count_byratio.txt
比较同源基因间这种调控作用是否存在差异
筛选出四组同源基因都有full-read 支持的基因对(表达了的同源基因对)
Ka/Ks :受到UTR调控的基因是否在序列进化水平上受到的选择压不一样。
python AddAnnotionTag.py all_gene_UTR_regulator.txt A2_D5_At_Dt_collinearity.txt homolog_UTR_regulator.txt
UTR调控在多倍化后丢失的
使用KOBAS3.0进行loacl GO 富集分析
提取基因的CDS序列,与拟南芥进行比对
筛选得分最高的拟南芥同源基因
使用拟南芥作为所有的基因集合作为背景基因集合,进行GO富集分析
##提取基因氨基酸序列,统一使用A2基因组的氨基酸序列
python ~/github/zpliuCode/Isoseq3/06UTRregulation/extractSequenceById.py /public/home/zpliu/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.pep.v1.0.fasta 01all_low_UTR.txt 01all_low_UTR.fa
##参考之前写的local KOBAS进行GO富集分析
与参考基因组编码框比较对转录本进行分类
in-frame change
frameshift (缺缺补补的不是3的倍数)
unchange ORF
other
for i in 1
do
##other类型
grep noreference ORF.txt2|wc -l
awk '$4==""{print $0}' ORF.txt2 |wc -l
##没有改变编码框的前提下,改变了氨基酸序列
grep inframeChange ORF.txt2|wc -l
awk '$9=="noframeshift"&&($4!=$7||$5!=$8){print $0}' ORF.txt2 |wc -l
##ORF没有改变
awk '$9=="noframeshift"&&$4==$7&&$5==$8{print $0}' ORF.txt2 |wc -l
##改变了编码框
awk '$9=="frameshift"{print $0}' ORF.txt2 |wc -l
done
基因组
inframeChange
frameshift
unchange ORF
Other
A2
20570
10975
23402
12165
D5
16877
7892
22578
4616
TM1
27011
13468(16.2%)
30082
12831(15.4%)
基因组
frame-shift
in-frame
Other
TM-1
13468(16.2%)
57093(68.4%)
12831(15.4%)
D5
7892(15.2%)
39455(75.9%)
4616(8.9%)
A2
10975(16.4)
43972(65.5%)
12165(18.1%)
针对ORF的分类去对Annotion与unAnnotion两种转录本进行注释
可以看出Annotion转录本的ORF大部分与参考基因组是一致的;而unAnnotion转录本的ORF大部分发生了frameshift和inframeChange
python AddAnnotionTag.py ../collapse/PacBio_Annotion/all_isoform.txt ORF.txt 11
for i in 1
do
##没有改变编码框的前提下,改变了氨基酸序列
awk '$NF=="Annotion"{print $0}' ORF_Annotion.txt|grep inframeChange |wc -l
awk '$NF=="Annotion"{print $0}' ORF_Annotion.txt|awk '$9=="noframeshift"&&($4!=$7||$5!=$8){print $0}' |wc -l
##ORF没有改变
awk '$NF=="Annotion"{print $0}' ORF_Annotion.txt|awk '$9=="noframeshift"&&$4==$7&&$5==$8{print $0}' |wc -l
##改变了编码框
awk '$NF=="Annotion"{print $0}' ORF_Annotion.txt|awk '$9=="frameshift"{print $0}' |wc -l
done
转录本类型
inframeChange
frameshift
unchange ORF
TM1_Annotion
13070
4933
25816
TM1_unAnnotion
13941
8535
4266
A2_Annotion
7106
3221
15484
A2_unAnnotion
13464
7754
7918
D5_Annotion
9805
3884
19872
D5_unAnnotion
7073
4008
2706
Annotion与unAnnotion两种转录本在编码框是否发生偏移上的卡方测验
转录本类型
frame-shift
other
TM1_Annotion
4933(10.7%)
41194
TM1_unAnnotion
8535(22.9%)
28730
A2_Annotion
3221(11.3%)
25323
A2_unAnnotion
7754(20.1%)
30816
D5_Annotion
3884(11.2%)
30945
D5_unAnnotion
4008(23.4%)
13133
两种转录本在终止密码子上的分析
由于选择不同的起始位点会对下游的终止密码子造成误判。
在分析终止密码子是否提前的时候,看那些起始密码子与参考转录本一致但是,终止密码子是否发生了提前。
从所有的转录本里一共筛选到
A2: 39947
D5: 35763
TM-1: 51052
##终止密码子没有改变
awk '$5==$8&&$6==$9{print $0}' 111 |wc -l
##终止密码子提前
awk '$1=="+"&&$6<$9{print $0}$1=="-"&&$5>$8{print $0}' 111 |wc -l
##终止密码子滞后
awk '$1=="+"&&$6>$9&&$5==$8{print $0}$1=="-"&&$5<$8{print $0}' 111 |wc -l
在这些转录本中,终止密码子的状态
TM1 17318 20.7%的转录本存在提取或者滞后的终止密码子
A2 14004 20.8%
D5 10290 19.8%
基因组
early
later
没有改变
TM1
7659
9659
33734
A2
7277
6727
25943
D5
5003
5287
25473
转录本类型
early+Later
没有改变
TM_Annotion
7458(21.4%)
27394(78.6%)
TM_unAnnotion
9860(60.9%)
6340(39.1%)
A2_Annotion
4681(23%)
15677(77.0%)
A2_unAnnotion
9323(47.6%)
10266(52.4%)
D5_Annotion
5412(20.2%)
21425(79.8%)
D5_unAnnotion
4878(54.6%)
4048(45.4%)
卡方测验:在unAnnotion转录本存在更高比例的提前终止的密码子
在多个转录本的基因中,转录本编码框没有改变的基因相比于转录本编码框改变的基因有更高的表达水平。
异常的剪接可能会调控基因的表达量的下调,对同源基因的剪接进行分析
同源基因间剪接的差异导致其表达量的差异。
##只转录出与参考转录本一致ORF的基因
awk '$4==$7&&$5==$8{print $0}' ../ORF_Annotion.txt |cut -f1|sort |uniq -c >ORF_same_with_reference.txt
cat ~/work/Alternative/result/Gh_result/CO31_32_result/evolution4/TM1_PacBio.txt |cut -f1|sort|uniq -c|cat - ORF_same_with_reference.txt |sort |uniq -d >11
##转录出与参考转录本一致的转录本同时也转录出提前终止的密码子转录本的基因
两种转录本表达量存在差异
从ORF的数据来看,unAnnotion转录本大部分发生了ORF的改变;同时两者的表达水平存在差异;这个分析可以说明编码框的改变可能会导致基因表达水平的改变;在At、Dt同源基因中
##提取转录本的表达水平
awk 'NR>=2{print $6"\t"$6"\t"$0}' /public/home/zpliu/work/Alternative/result/homologo/FEST3/geneExpress/stringtie/TM1/t_data.ctab >all_isoform_FPKM
##给转录本打上Annotion或者是unAnnotion的标签
python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py all_isoform.txt all_isoform_FPKM 11
mv 11 all_isoform_FPKM
##给所有转录本打上是否AS的标签
python ../ORF/AddAnnotionTag.py Alternative_isoform.txt ../collapse/PacBio_Annotion/all_isoform_FPKM all_isoform_FPKM
AS可能改变编码框使得转录本水平发生改变
##分析每种剪切中有多少比例的in-frame change,多少比例的early stop
python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py ../ORF/ORF.txt A3_isoform.txt A3_ORF.txt
##发生early stop的比例
for i in RI SE A3 A5
do
awk '$13=="earlyStop"{a+=1}END{print a/NR}' ${i}_ORF.txt
##发生in-frame change的比例
awk '$12=="noframeshift"&&($4!=$7||$5!=$8){a+=1}$12=="inframeChange"{a+=1}END{print a/NR}' ${i}_ORF.txt
done
找例子,比如一个基因存在很多 unAnnotion转录本,并且这个基因的的表达水平非常的低,同时也存在正常的转录本;就像nature plant那篇文章那样
做附表
awk '$9=="frameshift"{print $0}' ORF.txt2 |awk '{print "G. hirsutum\t"substr($1,1,8)"\t"$4"\t"$5"\t"$2"\t"$6}' >>~/work/Alternative/result/Gh_result/
两种转录本中AS事件的比例
基因组
Annototed
unAnnotated
P-value
A2
8146/28544
19865/38570
D5
17336/34829
9000/17141
TM1
16958/46127
20851/37265
Annotated isoforms 均值:0.285, 0.498, 0.368 平均 0.384
unAnnotated isoform 均值:0.515, 0.525, 0.560 平均0.533
找unAnnotion 基因进行基因的GO富集分析:
unAnnotion 转录本,同时导致提前终止的密码子,基因进行GO富集,
存在unAnnotion转录本的这些基因主要富集在一些与代谢和响应外界信号的一些通路中。
找到已经报导的基因的转录本模式图
##找基因
awk '$1=="+"&&$6<$9{print $0}$1=="-"&&$5>$8{print $0}' 111 |grep unAnnotion|cut -f2|sort|uniq |xargs -I {} grep {} 111 | awk '{print $2"\t"$NF}' |awk '$2=="unAnnotion"{a[$1][2]+=1}$2=="Annotion"{a[$1][1]+=1}END{for(i in a){print i"\t"a[i][1]"\t"a[i][2]}}'|awk '$2!=""&&$3!=""{print $0}'|less
比较同源基因的表达偏好性是否与unAnnotion转录本的数目存在相关性
Last updated
Was this helpful?