# 多倍化过程中的变化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
##比较亚组间转录本的保守性
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%的同源基因对**不存在功能保守的转录本**。

```bash
##都表达的同源基因中存在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值差不多。

```bash
##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

```bash
 ##获取特异性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甲基化数据的时候

```bash
##对筛选的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

```bash
 ##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
  >
  >   ```bash
  >   #根据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数目上的差异程度。

```bash
##保守的亚基因组同源基因，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的差异

```bash
##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`&#x20;
> > * `A2 D5 0.413`
> > * `At Dt 0.248`
>
> **A、D亚基因组；转录本功能上更加分化，但是AS剪接更加保守**。
>
> AS的保守与转录本的功能保守性之间是否存在联系（柱状图）
>
> **多少**比例的分化是由于AS导致的
>
> 讨论：这种保守的**AS信号可能是由于同处于一个细胞核的剪接复合体的识别**造成，三维结构上的靠近。
>
> * ~~存在保守AS的基因，转录本的保守性，与不存在保守AS的基因，转录本功能的保守性；**AS的保守性和转录本功能的保守性存在关联**~~
> * **转录本功能分化程度和基因AS剪接保守程度是否存在弱相关性**（柱状图）