多倍化过程中的变化3

分析同源基因对间转录本的保守性

统一用21066对四元同源基因对

在21066对同源基因中，转录出转录本的基因和没有转录本的基因数

  • A2 14602基因存在PacBio转录本

  • D5 13122基因存在PacBio转录本

  • AT 11543基因存在PacBio转录本

  • Dt 11749基因存在PacBio转录本

两个基因组间同时存在转录本的基因数

  • A2 vs D5 11885

  • A2 vs At 10670

  • D5 vs Dt 10198

  • At vs Dt 9063

四个同源基因都存在PacBio转录本：

8087个四组同源基因同时存在PacBio转录本

  • A2_vs_At

  • A2_vs_D5

  • At_vs_Dt

  • D5_vs_Dt

##比较亚组间转录本的保守性
bash Gh_result/05polyploidization/01isoformConserved/conserve_isoform.sh
##筛选保守的转录本
awk '$3~/A2/&&$4~/At/&&$6>=0.95{print $0}' Isoform_CDS.txt >A2_vs_At_conservedIsoform.txt

比较

存在保守转录本

不存在

A2 vs At

6112(57.3%)

4558(42.7%)

D5 vs Dt

6079(59.6%)

4119(40.4%)

A2 vs D5

7365(61.9%)

4520(38.1%)

At vs Dt

3249(35.8%)

5814(64.2%)

首先我们对同源基因都存在PacBio转录本进行了筛选，一共筛选得到了10670对A2 vs At和10198对D5 vs Dt同源基因对，

分析了各个同源基因内AS基因数目和AS事件：

  • A2和At同源基因分别有

    • A2 5430个基因存在16791个AS事件

    • At 4912个基因存在13628个AS事件

  • D5和Dt同源基因分别有

    • D5 5843个基因存在17136个AS事件

    • Dt 4745个基因存在13238个AS事件

比较了这些同源基因对转录本的保守情况，一共鉴定到

6112(57.3%)和6079(59.6%)存在保守转录本的同源基因对；同时我们发现30.4%-32.7%的同源基因对不存在功能保守的转录本。

##都表达的同源基因中存在AS的数目
A2_At_expression_A2_AS.txt
A2_At_expression_At_AS.txt
##

21066对同源基因中表达的基因数和AS事件数目

基因组

AS事件

基因数

A2

18814

6417

D5

18991

6805

At

13981

5122

Dt

14083

5115

多倍化过程中保守的AS比例：

(2885+3418)/(18814+18991)

分析同源基因间保守的AS事件：

使用K-mer获得定位另外一个基因组上的坐标，然后将k-mer坐标与AS坐标取交集，获取保守的剪接事件。

  • k-mer坐标至少需要90%与AS坐标存在交集

  • 所有的AS事件都保守

  • 存在特异性的AS事件

分析保守的剪接事件之间序列的相似程度是否比较高

特异性的剪接事件对应的PSI值差不多。

##D5 vs Dt
 cat splitFile/*kmer.txt >AS_kmer.txt
 ##获取k-mer对应的剪接事件
 python ~/github/zpliuCode/Isoseq3/04ASconserved/extractConservedAS.py  AS_kmer.txt  ../../../05polyploidization/02ASevent/Dt_homolog_AS.txt 11
 awk '$2!~/>/{print $0}' 11  >D5_Dt_conservedAS.txt

获取特异性的剪接事件对应是constitutive Exon或者Intron

 ##获取特异性IR和SE事件对应k-mer为constitutive-exon或者intron
 awk '$2~/>/{print $0}' 11|grep RI >D5_specialRI.txt
 awk '$2~/>/{print $0}' 11|grep SE >D5_specialSE.txt
 ##提constitutive exon
 intersectBed -loj -a gene.bed AS_length/constitutive_exon.bed -F 1|awk '{OFS="\t";print $7,$8+1,$9-1,$4,$6,"constitutiveExon"}' >TM1/constitutiveExon.bed
 ##提取constitutive intron
  intersectBed -loj -a gene.bed AS_length/constitutive_intron.bed -F 1|awk '{OFS="\t";print $7,$8,$9,$4,$6,"constitutiveExon"}' >TM1/constitutiveExon.bed
  ##分类
  python ~/github/zpliuCode/Isoseq3/04ASconserved/extractSpecialAS.py D5_specialRI.txt ../../constitutiveFeature/TM1/constitutive_feature.bed  3
  ##计算每种事件的PSI值
  python ~/github/zpliuCode/Isoseq3/03ASquantify/extractASPSI.py all_AS_quantify.txt ASeventFile.txt outFile.txt

分析存在保守AS事件的同源基因

同源基因间存在保守的AS数和基因数,

A基因组间存在保守AS的基因比例比D基因组小，A亚组进化快一些

AS的保守程度：（那种一对多的情况，如果大于1就当做1来算）

都存在AS的同源基因对，比较这些同源基因的AS保守情况

(保守的AS数a+b)/(同源基因的AS数 a+b)

  • A基因组和D基因组在多倍化前后

  • A基因组AS更保守还是D基因组

D基因组AS相比于A基因组更保守，多倍化后At、Dt间AS剪接更加保守

RI和SE在多倍化过程中的变化

A2 vs At

保守的RI事件：

保守的SE事件：

A基因组中特异性AS事件的去向

type

RI

SE

exon

intron

Other

Special Total

A2RI

1481

-

530 (7.6%)

3984 (57.2%)

2457(35.3%)

6971

A2SE

-

182

295(33.8%)

325(37.3%)

251(28.9%)

871

AtRI

1481

-

439(9.6%)

2379(52.0%)

1758(38.4%)

4576

AtSE

-

182

399(42.0%)

207(21.8%)

345(36.2%)

951

D基因组中特异性AS的去向

type

RI

SE

exon

intron

Other

special Total

D5RI

1410

-

576(9.7%)

3070(51.4%)

2322(38.9%)

5968

D5SE

-

287

873(55.6%)

319(20.3%)

378(24.1)

1570

DtRI

1410

-

523(11.0%)

2380(50.2%)

1839(38.8%)

4742

DtSE

-

287

284(35.9%)

258(32.6%)

249(31.5%)

791

多倍化过程中大量亚组特异性的AS事件，转变为了组成型的intron和组成型的exon；

比较多倍化前后保守AS的PSI和PSJ变化

由于测序read覆盖不足的原因，鉴定出的很多AS可能是短暂出现，可能是剪接机制一不小心错误的识别了冒名顶替的剪接位点。为了准确的鉴定多倍化过程中AS的变化，我们对物种中包含多转录本的基因的AS event进行了量化：

过滤掉那种只由参考转录本和PacBio比较得出的AS，这种AS无法进行量化，对于与参考转录本比较得到的AS时，参考转录本的reads数默认为0

$a/(a+b+c)$

a: 包含剪接事件的转录本对应的FL read数目

b、c: 不包含剪接事件的转录本对应的FL read 数目

junction reads that support inclusion of alternative

junction reads that support exclusion of alternative exons

对于RI事件，支持的read数目还要进行标准化

read count read length 2/(read length *2+event length)

  • PSI 包含这个剪接事件的比例

    对于AS变成constitutive exon或者constitutive intron；这个也许由于测序深度的影响造成AS没有检测到，分析PSI值差异大不大。

  • constitutive exon PSI为1

  • constitutive intron PSI为0

    PSI再乘以基因的表达量，之后用于DNA甲基化数据的时候

##对筛选的AS事件，进行量化03quantifyASevent
python ~/github/zpliuCode/Isoseq3/03ASquantify/caculateAS_ratio.py  ../02ASevent/A2_homolog_AS.txt  ~/work/Alternative/result/Gh_result/CO31_32_result/evolution4/A2_PacBio.txt  11
##使用RNA-seq的数据进行量化
python ~/github/zpliuCode/Isoseq3/03ASquantify/PSIByRNAseq.py  1  ~/work/Alternative/result/homologo/FEST3/geneExpress/hisat2/D5_rmdup.bam 22

由多倍化过程中内含子转变为IR的PSI值明显的低于外显子转变为IR，存在剪接印记效应，这种剪接转录本的含量比较低，或者是发挥完作用后被降解的过程中被测到了。

为了找到多倍化过程中稳定且保守的RI事件。我们使用另外一个RNA-seq的短读read数据对RI事件的含量进行了量化。同样发现了类型的现象，根据short-read得到稳定的RI事件。

  • IR2intron

  • IR2exon

  • intron2IR

  • exon2IR

  • SE2intron

  • SE2exon

  • exon2SE

  • intron2SE

  • IR2IR

  • SE2SE

IR转变为intron的例子中，IR的PIR值很低，这会不会是事件只是短暂的出现了一会

虽然这个地方检测到了RI的存在，但主要还是以非RI的转录本为主。

但是也存在打破这种诅咒的现象；用RNA-seq做一个例子的图

  • intron变成IR后，IR的PSI值很高，read数目也很多

分析RI和SE在四个亚基因组中的保守性

保守的RI和SE事件

  • A2 D5 1418 180

  • At Dt 1107 179

  • A2 At 1491 182

  • D5 Dt 1412 287

 ##AS在四个基因组中的保守性
 python ~/github/zpliuCode/Isoseq3/04ASconserved/homologConservedAS.py  all_conserve_RI_SE.txt  11
 ##筛选同一个位置处得到最保守的剪接事件，
  python ~/github/zpliuCode/Isoseq3/04ASconserved/filterConservedAS.py  ../A2_D5/A2_D5_conserved_AS.txt 11
 ##二倍体和四倍体中保守的AS事件
 awk '$1~/RI/||$1~/SE/{print $0}' ../A2_D5/A2_D5_conserved_AS.txt |wc -l
 ##A基因组和D基因组保守的AS事件

  • 在二倍体和四倍体都保守的IR和SE事件

都保守的AS事件：577

只在二倍体中存在，2274

只在四倍体中存在 1981

  • 在A基因组和D基因组保守的IR和SE事件

只在A基因组中存在，2308

只在D基因组中存在 ：2841

亚组特异性的剪接事件？

举例子:

总结：

  • 分析存在保守AS的基因占（存在AS基因）的比例（柱状图）

  • 分析多倍化过程中IR向intron和exon类型转变（很多特异AS事件转变为CI）

  • 不同转变类型的PIR值大小（箱型图），特例图（打破诅咒的）

  • 比较不同基因组间AS的保守程度（箱型图）

  • 二倍体和四倍体都存在的保守IR和SE（韦恩图,胶图）

    举个例子（PacBio转录本例子）

    筛选标准：

    • read数目都大于10

    • PSI值都大于0.1

    • 保守的RI事件之间，PSI值相差两倍以上

    • 不保守的RI事件，则是PSI值介于0.05-0.95，支持read>5

      #根据PSI值筛选差异的RI事件，进行GO富集分析
paste A2RI_PSI.txt AtRI_PSI.txt |awk '$7-$3>$3||$3-$7>$7{print $0}'|less 
##特异性的RI事件
awk '$3>=0.05&&$3<=0.95&&$2>5{print $0}' A2RI_PSI.txt |less

    在这些保守的IR或者SE事件中,有多少PSI值是存在两倍差异的

基因组

二倍体上调

四倍体上调

不变

A RI

1481

96()

553()

832()

D RI

1410

106()

345()

959()

A SE

25()

21()

136()

D SE

287

35()

21()

231()

对差异表达的保守RI事件进行GO富集分析

第五个结果：多倍化过程中转录本的保守性

模型：AS加速同源基因间的进化

比较转录本的趋同进化和并行分化

2418个同源基因在多倍化前后都存在功能保守的转录本。

多倍化前有多少比例的同源基因存在保守转录本：

  • 二倍体中 7365/11885 (62.0%) 基因存在保守转录本

  • 四倍体中 有3249/9063 (35.8%)基因存在保守转录本

  • A2 vs At中有6112/10670 (57.3%)基因存在保守的转录本

  • D5vsDt中有6079/10198 (59.6%) 基因存在保守的转录本

比较不同基因组间转录本的保守程度

  • 多倍化前后，不同基因组存在保守转录本的基因的比例变化（柱状图）

  • 分析四个同源基因间存在保守转录本的基因数（韦恩图）

  • 分析存在保守转录本的同源基因，相比于功能发生分化的同源基因表达水平和Ka/Ks值（箱型图）

多倍化过程中存在保守转录本基因与非保守基因表达水平和Ka\Ks值

并行分化

虽然At、Dt间保守转录本的程度比较低，只有（3249）30%多。5814个亚基因组同源基因功能发生了分化，在这些功能分化的同源基因中有80.3%的基因是与二倍体祖先存在功能一致的转录本的。

  • 2346 个At基因与A2存在保守的转录本 (饼状图)

  • 240个At基因与D5存在保守的转录本

  • 2322个Dt基因与D5存在保守的转录本

  • 231个Dt基因与A2存在保守的转录本

在At与Dt转录本分化的这些同源基因中，有基因是分别转录了二倍体基因的不同转录本，基因的AS直接变成不同基因来转录不同的转录本，形成功能互补。总的来说一共有4402基因与二倍体相比存在保守的转录本，而亚基因组同源基因间不存在保守的转录本。

这些同源基因转录本发生分化的基因中，基因表达水平和Ka/Ks值的差异（柱状图）

比较分化的At、Dt基因与保守的转录本的At、Dt基因，在AS数目上的差异程度。

##保守的亚基因组同源基因，AS数目上的差异程度
awk 'NR%2!=0{printf $0"\t"}NR%2==0{print $0}'  conserved_AScount.txt |awk '{print sqrt(($2-$4)*($2-$4))}' >1
##不保守的亚基因组同源基因，AS数目上的差异程度
awk 'NR%2!=0{printf $0"\t"}NR%2==0{print $0}' diversity_AScount.txt |awk '{print sqrt(($2-$4)*($2-$4))}' >2

## 计算相差AS数目对应的频率
sort 1 -n|uniq -c |awk '{print $1/3249"\t"$2}'

GO富集图（分析）

举例子（igv画转录本结构图）

在发生分化的At、Dt同源基因中，有多少对同源基因是祖先同源基因功能分化的结果，两个同源基因分别转录出祖先不同的转录本，承担不同的功能，亦或者是通过AS来调节同源基因的表达。

  • 在2346 个At与A2存在保守转录本的基因中，有296个Dt基因与A2相比也存在保守的转录本，

  • 在2322个Dt与D5存在保守转录本的基因中，有221个At基因与D5相比存在保守的转录本。

有的At、Dt同源基因既可能继承了A2的互补转录本，也可能继承了D5的互补转录本；所以去个重，总共得到了464个同源基因对，它们分别转录出祖先基因中的一个转录本，功能上是互补。

  • 同源基因产生功能互补的转录本，这些转录本存在AS的差异

##At、Dt亚组间转录本功能发生分化的基因，但是都与祖先基因组存在保守的转录本 ：/public/home/zpliu/work/Alternative/result/Gh_result/05polyploidization/01isoformConserved/homologDifferentiation
cut -f3 At_Dt_differ_D5_Dt_gene.txt |xargs -I {} grep {} ../D5_vs_At/conserved_Isoform.txt |cut -f1 |sort |uniq |xargs  -I {} grep {} ~/work/Alternative/result/Gh_result/CO31_32_result/evolution5/A2_D5_At_Dt_collinearity.txt  >At_DtDifferent_AllConservedWithDiploid/At_Dt_different_AllConservedWithD5.txt
##对应的转录本，在A2中是否存在AS差异
python ~/github/zpliuCode/Isoseq3/08IsoformDifferentiation/homologIsoformDifferentiation.py ~/work/Alternative/result/Gh_result/CO31_32_result/evolution2/D5_AS.txt  ../../D5_vs_Dt/D5_vs_Dt_conservedIsoform.txt ../../D5_vs_At/conserved_Isoform.txt 22 At_Dt_different_AllConservedWithD5.txt

转录本的保守程度（箱型图）

A基因组的转录组经历更大变化

不同基因组间转录本的保守程度:

保守转录本数目/(转录出转录本的数目) 保守程度：

  • A2 At 0.386

  • D5 Dt 0.421

  • A2 D5 0.413

  • At Dt 0.248

A、D亚基因组；转录本功能上更加分化，但是AS剪接更加保守。

AS的保守与转录本的功能保守性之间是否存在联系（柱状图）

多少比例的分化是由于AS导致的

讨论：这种保守的AS信号可能是由于同处于一个细胞核的剪接复合体的识别造成，三维结构上的靠近。

  • 存在保守AS的基因，转录本的保守性，与不存在保守AS的基因，转录本功能的保守性；AS的保守性和转录本功能的保守性存在关联

  • 转录本功能分化程度和基因AS剪接保守程度是否存在弱相关性（柱状图）

Previous设计引物Next不同棉种间AS的差异

Last updated