多倍化过程中的变化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剪接保守程度是否存在弱相关性（柱状图）