可变剪切的进化分析

对同源基因的AS和ORF进行统计

对于AS isoform的分类不能用AS来描述了；改成使用是否包含这段AS片段

例如Ghir_A01G005170基因中对应的PB.5574.2转录本就存在一个ES事件,但是这个转录本不包含对应的exon序列；suppa软件的输出中将其放在最后一列。

对每个棉种中基因的AS和ORF数据进行统计分析：

• 基因存在发生AS的转录本，与正常功能的转录本；并且对应的ORF发生了提前终止、或者frameshift

• 基因存在发生AS的转录本，与正常功能的转录本；并且对应的ORF编码框没有发生改变，发生了in-frame change；在中间增加了一个序列；或者later stop

梳理一下AS的结果

A3:Ghir_A01:116840-116977:116835-116977:-

两个转录本上的exon坐标前面表示比较长的那一段exon序列；后面那一段表示比较短的那段exon序列

RI:Ghir_A01:524019:524047-524087:524141:-

表示两端exon的坐标序列；从第二个到第三个则是对应保留的intron坐标，第一个转录本表示包含保留的Intron

SE:Ghir_A01:2306968-2307075:2307164-2307243:+

中间的坐标表示被跳跃的exon坐标；两端的坐标分别表示两端exon的位置，第一个转录本表示包含跳跃的exon

A5:Ghir_A01:242983-244582:242983-244631:-

前面的坐标表示包含A5那段序列的转录本，在两端exon的坐标；后面则是不包含A5那段序列的转录本

针对ORF和AS进行分类

• 转录本与注释的转录本再ORF上存在差异、转录本与对应的注释转录本再AS上同样存在差异

• 转录本与注释的转录本ORF上不存在差异，转录本与对应的注释转录本再AS上不存在差异

• 转录本与注释的转录本ORF上不存在差异，转录本与对应的注释转录本再AS上存在差异

• 转录本与对应的注释转录本存在差异，但是与对应的注释转录本再AS上不存在差异

对转录本进行打上标签：

• 相比于对应的参考转录本，发生了AS，并且包含剪切片段

• 相比于对应的参考转录本，发生了AS，并且不包含剪切片段

• 相比于对应的参考转录本，没有发生AS

## 提取所有AS的信息
awk 'NR>=2{print $0}' D52_A5_strict.ioe  >>All_alternative.txt 
##给PacBio转录本进行注释
python  ~/work/Alternative/result/Gh_result/CO31_32_result/AS2/Add_AsTage.py All_alternative.txt  ../ORF/ORF.txt  11

根据是否与AS相关联将PacBio转录本分成两类：

比较noAS 转录本与AS转录本在移码框和终止密码子上是否存在差别；然后再看每种剪切事件对应的转录本有多少比例发生移码突变，有多少发生终止密码子提前

• 分析与AS相关的转录本的组成

  将isoform associate with Alternative splice 转录本分为两类：

  1. contain alternative splice fragment 相比于参考基因组转录本

  2. not contain alternative splice fragment

基因组	与AS相关	与AS无关
TM1	24407	46154
A2	13547	41398
D5	13433	33915

统计AS associate 相关的isoform中多少AS events包含了alternative splice fragment

可以发现在大多数IR事件中，PacBio检测到的AS associate isoform是包含了对应的intron的；而A3、A5、ES事件中；包含AS片段的转录本比例与不包含的比例接近1:1；这也可能是导致PacBio转录本相比于参考基因组的cDNA序列长度更长的原因之一；

基因组	RI	ES	A3	A5
TM1	9069/1673	1231/1371	6338/4507	4822/3752
A2	6186/443	726/568	2970/2097	2380/1674
D5	5421/630	785/673	2767/2807	1821/1561

• 分析AS associate isoform与编码框的偏移

for i in 1
do
grep  NoAS 11 |awk '$9=="frameshift"{print $0}'|wc -l
grep  NoAS 11 |awk '$9=="noframeshift"{print $0}'|wc -l
grep  NoAS 11 |awk '$9=="inframeChange"{print $0}'|wc -l
done
##移码框的卡方检验
>c=matrix(c(5583,7885,40571,16522),nrow=2)
>chisq.test(c) ##卡方检验
>fisher.test(c) ##使用fisher T检验

基因组	frameshift	noframeshift	inframeChange
TM1_noAS	5583 12.10%	39621 85.85%	950 2.06%
TM1_AS	7885 32.31%	13820 56.62%	2702 11.07%
A2_noAS	5994	34402	1004
A2_AS	4981	7029	1537
D5_noAS	3981	29100	834
D5_AS	3911	7461	2061

• 分析AS associate isoform 与终止密码子的提前

##将AS associate isoform分类
for i in 1
do
grep -v NoAS 11 |awk '$10=="earlyqStop"{print $0}'|wc -l
grep -v NoAS 11 |awk '$11=="laterStop"{print $0}'|wc -l
grep -v NoAS 11 |awk '$9=="noframeshift"&&$5==$8&&$4==$7{print $0}'|wc -l
grep -v NoAS 11 |awk '$9=="inframeChange"{print $0}'|wc -l
grep -v NoAS 11 |awk '$5==$8&&$4!=$7{print $0}'|wc -l
done
## 卡方检验终止密码子是否发生改变
>c=matrix(c(8444,9910,37710,14497),nrow=2)
>chisq.test(c) ##卡方检验

基因组	early	later	both unchange	ORF change but stop not change
TM1_noAS	7420	1024	27669	9091+950
TM1_AS	8523	1387	2413	9382+2702
A2_noAS	8476	425	23395	8100+1004
A2_AS	6203	67	7	5733+1537
D5_noAS	5713	1009	21102	834+5257
D5_AS	4763	519	1476	2061+4614

分析四种剪切事件中发生frameshift和early stop codan的比例

##只分析那种存在一种剪切事件的转录本
grep -v NoAS 11 |awk '$9=="inframeChange"&&$11==""{print $0}$9!="inframeChange"&&$15==""{print $0}'
##A5发生frameshift的数目
for i in 1
do
grep -v NoAS 11 |awk '$9=="inframeChange"&&$11==""{print $0}$9!="inframeChange"&&$15==""{print $0}'|grep ":RI:" |awk '$9=="frameshift"{print $0}'|wc -l  
##A5总数
grep -v NoAS 11 |awk '$9=="inframeChange"&&$11==""{print $0}$9!="inframeChange"&&$15==""{print $0}' |grep ":RI:"|wc -l
done
##A5发生early stop的数目
grep -v NoAS 11 |awk '$9=="inframeChange"&&$11==""{print $0}$9!="inframeChange"&&$15==""{print $0}'|grep ":RI:" |awk '$10=="earlyStop"{print $0}'|wc -l

进行卡方检验

编码框的改变

类别	发生frameshift	不发生frameshift	p-value
TM1_noAS	5583	40571	2.2e-16
TM1_AS	7885	16522
A2_noAS	5994	35406	2.2e-16
A2_AS	4981	8566
D5_AS	3981	29934	2.2e-16
D5_noAS	3911	9522

终止密码子的改变

类别	终止密码子发生改变	终止密码子没发生改变	p-value
TM1_noAS	8444	37710	2.2e-16
TM1_AS	9910	14497
A2_noAS	8901	32499	2.2e-16
A2_AS	6270	7277
D5_noAS	6722	27193	2.2e-16
D5_AS	5282	8151

发生frameshift的比例

AS couple with frame shift and early stop

ES事件导致的frameshift比例是最低的可能和exon是3的倍数相比于intron更高有关；这里可以统计一下

基因组	IR	ES	A3	A5
TM1	1822/6048	401/1492	1734/6157	1327/4387
A2	1428/4234	298/850	1159/3148	937/2457
D5	1187/3952	223/978	938/3699	615/2249
average	0.312	0.279	0.300	0.318

发生early stop的比例

基因组	IR	ES	A3	A5
TM1	2577/6048	350/1492	1588/6157	1267/4387
A2	2293/4234	338/850	1184/3148	975/2457
D5	1899/3952	232/978	933/3699	552/2249
average	0.483	0.29	0.295	0.31

小结

在证明了编码框的改变会影响基因表达后；我们对可变剪切对ORF的影响进行了探究；首先将PacBio与参考转录本进行比较鉴定出isoforms associate with Alternative splice；和none Alternative splice；其中IAAS转录本又分为两组：

• 与参考转录本相比，包含剪切片段的

• 与参考转录本相比，不包含剪切片段的

通过统计这两种AS 转录本在各种剪切事件的比例；发现在IR事件中包含剪切事件片段的转录本占据主要；而其他事件中两种转录本的比例约为1:1；

进一步比较noneAS转录本与AS associate 转录本在frameshift、earlyStop上是否存在差异；进行卡方测验发现AS associate 与非AS associate isofrom转录本在编码框的改变、提取终止的密码子上没有显著性差异

论证AS与编码框、提前终止的密码子存在联系后；提取AS转录本的子集，只看那些与参考基因组转录本存在一个剪切事件差异的剪切事件，分析这些剪切事件导致frameshift、earlyStop的比例；分析发现exon skip事件中导致framshift、earlystop的比例是最小的，而IR事件有着高比例的frameshift和earlyStop；这可能和exon原本就起着重要的作用

AS与基因表达

前面已经讨论了AS会改变转录本的编码框，并且ES事件对编码框的影响是最小的；而IR事件对编码框的影响是最大的；既然AS影响了转录本的编码框；在后续的翻译阶段就可能会存在着两种命运；

• 例如AS导致的编码框的移码、终止密码子的提取，最终受到NMD调控机制调节从而导致转录本被降解，这种机制在一定程度上相当于转录本的一个分流器，调控基因的表达

• AS保留或者截取的片段，并没有导致mRNA被降解的命运；而是在经过翻译后形成了新的蛋白质亚基，行使新的功能。

为此我们将经过AS和ORF注释后的转录本进行了分类：

• none AS isoform

• AS association isoform with (frameshift|earlyStop)

• AS association isoform with(inframeChange)

## none AS isoform
grep NoAS 11 >NoAS_ORF.txt
## 发生frameshift或者earlyStop的转录本
grep NoAS 11 -v |awk '$9=="frameshift"||$10=="earlyStop"{print $0}' >frameshift_earlyStop_ORF.txt
## 发生inframeChange的转录本
grep NoAS 11 -v |awk '$9=="inframeChange"{print $0}' >inframeChange_ORF.txt
##提取对应的表达量
python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py  ~/work/Alternative/result/homologo/FEST3/geneExpress/stringtie/D5/t_data.ctab  frameshift_earlyStop_ORF.txt frameshift_earlyStop_FPKM.txt

从整体表达水平看，None AS转录本的表达水平高于 inframeChange isoform转录本，而最小的则是frameshift转录本。从表达水平上的角度也验证了我们的想法。根据这个想法我们把从中筛选了基因的子集：

• 只转录出与参考转录本ORF一致转录本

• 转录出与参考转录本ORF一致的同时，还转录出一些由于AS导致frameshift|earlyStop的转录本

• 转录出与参考转录本ORF一致的同时，还转录出一些由于AS导致inframeChange的转录本；编码框没有变，ORF改变了

这里画一个韦恩图，多少基因转录出与参考转录本ORF一致的 转录本，多少基因转录出可能有AS导致的frameshift&earlyStop ，多少基因转录出in-frame Change；经过分析将近有18%的基因在转录出参考转录本的同时，还受到AS的调控；11.6%的基因在没有转录出于参考转录本一致的ORF情况下受到AS的影响。

##只提取转录出与参考转录本ORF一致的基因
awk '$9=="noframeshift"&&$5==$8&&$4==$7{print $0}' NoAS_ORF.txt |cut -f1|sort |uniq  >No_AS_genelist

for i in 1
do 
##转录出NoAS的同时受到AS调控
cat No_AS_genelist frameshift_genelist |sort |uniq -d >1
cat No_AS_genelist inframeChange_genelist |sort|uniq -d >>1
sort 1 |uniq >gene_refORF_AS
##只转录出NoAS的基因
cat No_AS_genelist gene_refORF_AS |sort |uniq -u >only_NoAS
##只受到AS调控
cat No_AS_genelist frameshift_genelist |sort |uniq -d|cat - frameshift_genelist |sort |uniq -u >only_AS
cat No_AS_genelist inframeChange_genelist |sort |uniq -d|cat - inframeChange_genelist |sort |uniq -u >>only_AS 
sort only_AS |uniq >1
mv 1  only_AS
rm 1 -rf
done

对基因分好类后比较它们的表达水平和进行GO富集分析

##提取它们的表达水平差异
 python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py  ~/work/Alternative/result/homologo/FEST3/geneExpress/stringtie/A2/gene_FPKM.txt only_NoAS only_NoAS_FPKM
##进行GO富集分析
python ~/work/Alternative/result/Gh_result/CO31_32_result/ORF/AddAnnotionTag.py ~/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_Noscagenes_GO_V3.annot only_AS only_AS.GO