可变剪切的翻译分析

分析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转录本的这些基因主要富集在一些与代谢和响应外界信号的一些通路中。

找到已经报导的基因的转录本模式图

• Ghir_D09G000560

  https://www.nature.com/articles/srep32196

  调控 PIN(生长素转运蛋白)的表达，从而控制生长素的分布与运输；调控拟南芥的生长发育

##找基因
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转录本的数目存在相关性