可变剪切的进化分析

对同源基因的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
Previous可变剪切统计Next保守AS模式的鉴定

Last updated