# 02

## 保守的 isforms鉴定

使用ab-Blast进行鉴定

#### 使用xdformat建库

* `-n`指定序列类型
* -o 输出数据库前称
* 最后跟fasta文件

```bash
xdformat  -n  -o test test.fa
```

#### 进行比对

参考说明 <https://blast.advbiocomp.com/doc/parameters.html>

* 库前称
* 序列文件
* `-o`输出文件

```bash
~/software/wu-blast/ab-blast-20200317-linux-x64/ab-blastn test test.fa -o 11  e=1e-5 mformat=2
```

#### 提取isform每个剪切位点附近的exon序列

* 提取exon两端锚点附近300bp序列
* 不足300bp则提取整个exon

**运行wu-blast**

* 同源基因都存在isform的情况

```bash
python ~/scripte/Alternative/module/FEST3/extract_spliceSiteSeq.py -p geneExpress/isform.gtf  -r gene_isformCount.txt2  -g ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta -o tmp
## 建库
~/software/wu-blast/ab-blast-20200317-linux-x64/xdformat -n -o all_vs_all all_vs_all.fa
##鉴定同源基因间保守的isform

python ~/scripte/Alternative/module/FEST3/wu-Blast.py  -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt  -r1 ../A2/gene_isformCount.txt2  -r2 ../TM1/At_isformCount.txt  -fa1 A2_isform_splice.fa  -fa2 At_isform_splice.fa  -o hhhhhh
##
```

#### 鉴定保守isform2.0

* exon数目相同，
* 从5‘到3’端每个exon的长度都接近

```bash
python ~/scripte/Alternative/module/FEST3/conserveIsform.py  -r1 ../A2/gene_isformCount.txt2  -r2 ../D5/gene_isformCount.txt2 -gtf1 ~/work/Alternative/result/Ga_result/CO11_12_result/07_annotation/merge.gtf -gtf2 ~/work/Alternative/result/Gr_result/CO41_42_result/07_annotation/merge.gtf -gap 15 -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_D5_collinearity.txt  -o zzz
```

### 鉴定保守的isform3.0

* wu-blast结果存在多对多的情况
* 比较isform的长度，取最接近的那个
* 去除isform长度差异不超过200bp的

```bash
##转换wu-blast结果
awk '$1~/evm/{print $0}$1~/Ghir/{print $2"\t"$1}' ../A2_At.txt
##提取保守的isform
python ~/scripte/Alternative/module/FEST3/conserveIsform2.py  -blast A2_At.txt  -r1 ~/work/Alternative/result/Ga_result/CO11_12_result/07_annotation/merge.gtf  -r2 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf  -o 222
##去重
sort 222 |uniq |sort -k2,2 -k5,5n|awk '{print $1,$3,$4,$5,$2}' OFS="\t" |uniq -f4|awk '{print $5,$2,$3,$4,$1}' OFS="\t" |sort -k5,5 -k4,4n|uniq -f4 |awk '{print $5,$1,$2,$3,$4}' OFS="\t" >A2_At.txt
## 筛选并且更改格式
awk '$5<=200{split($1,a,"-");split($2,b,"-");print a[1]"\t"a[2]"\t"b[1]"\t"b[2]}' ../conserve_isform/conserve_isform/D5_Dt.txt
```

#### 保守的isform统计

| 比较    | 都存在isform的基因数 | isform数 | 保守isform数 | 基因数  |
| ----- | ------------- | ------- | --------- | ---- |
| A2 D5 | 13543         | 77180   | 11823     | 8196 |
| At Dt | 10368         | 50282   | 6331      | 4850 |
| A2 At | 12041         | 68419   | 12868     | 8604 |
| D5 Dt | 12580         | 63143   | 12335     | 8570 |

统计存在保守isfrom的数目分布，以及对应的gene isform的总数的平均值

```bash
python ~/scripte/Alternative/module/FEST3/isformlevel/isformNumber.py -c conserve_isform.txt  -all isform_count.txt  -f1 0 -f2 2 -o1 test -o2 test2
```

![NN2lIP.png](https://s1.ax1x.com/2020/06/23/NN2lIP.png)

二倍体到四倍体，多倍化过程中，存在保守isform的基因，非保守的isform减少，表明isform的多态性在减少

```bash
python ~/scripte/Alternative/module/FEST3/isformConserveRate.py -all isform_count.txt  -c conserve_isform.txt  -o 11
```

#### 二倍体内保守，四倍体内不保守例子

isoform在多倍化后的丢失

* A2 与D5保守，只有At存在保守的
* A2与D5保守，只有Dt存在保守的

#### 二倍体内部保守，四倍体内保守的例子

isoform多倍化后，在两个亚基因组间共享

#### 二倍体内保守，四倍体内同样保守

```bash
python ~/scripte/Alternative/module/FEST3/isformlevel/polypolid.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_D5_At_Dt_collinearity.txt  -A2D5 ../A2_vs_D5/conserve_isform.txt  -AtDt ../At_vs_Dt/conserve_isform.txt  -A2At ../A2_vs_At/conserve_isform.txt  -D5Dt ../D5_vs_Dt/conserve_isform.txt  -o AtDt分化 -o1 AtDt趋同 -o2 保持不变
```

* 输出的结果还需要去一下重,由于在用wu-blast直接比对的时候没有鉴定出保守性，但是通过两个基因组间的关系间接证明了保守性，这个先去掉，原始文件先存着；到时候设计引物的时候在看看

```bash
## 多倍化后isoform丢失
awk '{print $0"\t"$4"\t"$5}' 111 |sort -k7,8 |uniq -f 6 -u |cut -f1-6 >polploid_decrease.txt
## 多倍化后趋同
 sort -k5,6 22 |uniq -f4 -u >polploid_inrease.txt
## 多倍化后仍旧保守的
mv 333  polploid_unchange.txt
```

## 分析四倍体类保守转录本的组成

At与Dt是保守的，但是只在A2中保守，或者只在D5中保守；趋同进化

\| At 与Dt 保守的6331 对isform中，有3166在二倍体内同样保守，1145对只在A2中出现，840对只在D5中出现，剩下1180对，其他

分析A2与D5中保守，但只在At与Dt一个亚基因组中保守的，并行分化

\| A2与D5中保守的11823条转录本，在四倍体中有3166条同样保守，有2262条在At中保守，2380条在Dt中保守，4,015是其他类型

**多倍化过程中两股力量同时推动着棉花的进化**

虽然At与Dt比较中存在较少比例的保守isofrom，但是我们分别对A2与D5中保守的isofrom，与两个亚基因组比较后发现39%，4600多条只在At或者只在Dt中保守的isofrom；与此同时通过对At与Dt中已经鉴定到的保守的isform，我们发现有将近31%的isofrom只在A2或者只在D5中保守，更多的D基因组向A基因组方向进化

![NqjSq1.png](https://s1.ax1x.com/2020/07/02/NqjSq1.png)

#### 二倍体A2、D5合并后的isform数目与At、Dt合并后是否存在差异

* At相比于A2 isform数减少
* Dt相比于D5 isform数减少
* At与Dt合并后相比于A2、D5合并后isform数减少
* A2 与D5 间存在差异而在At、Dt中没有差异

![NDvp59.png](https://s1.ax1x.com/2020/06/26/NDvp59.png)

#### 保守的isform表达水平与非保守isform表达水平比，

> A2 conserve vs A2 uniq | At conserve vs At uniq
>
> 保守的isoform表达水平更高

```bash
##获取每种比较中isform和对应的fpkm值
python ~/scripte/Alternative/module/FEST3/isformlevel/isformFPKM.py  -c conserve_isform.txt  -all isform_count.txt  -fpkm1 ../../geneExpress/stringtie/A2/t_data.ctab  -fpkm2 ../../geneExpress/stringtie/D5/t_data.ctab   -o A2_D5_conserveFPKM.txt  -uniq A2_D5_uniqFPKM.txt
```

![NyQLfH.png](https://s1.ax1x.com/2020/06/27/NyQLfH.png)

#### 四组同源基因中保守的isform表达水平的变化是否存在剂量效应

#### 保守的isform与非保守的isform发生AS的频率

```bash
##获取发生剪切的isform
awk '{print $($8)}' ~/work/Alternative/result/Gr_result/CO41_42_result/11_AS/end_splice.txt3 |sort |uniq >splice_isform.txt
##提取每种类型中剪切的数目
awk '$1~/Ghir/{print $2}' A2_At_uniq_FPKM.txt |xargs  -I {} grep -E "{}$" ../TM1/splice_isform.txt |wc -l
```

| 比较       | 保守isform    | 发生剪切isform | 不保守isform   | 发生剪切isform |
| -------- | ----------- | ---------- | ----------- | ---------- |
| A2 vs At | 12868/12868 | 1445/1690  | 23680/13270 | 7232/4680  |
| D5 vs Dt | 12335/12335 | 990/1471   | 17868/14053 | 5305/4817  |
| A2 vs D5 | 11823/11823 | 1121/841   | 21264/15536 | 6533/4750  |
| At vs Dt | 6331/6331   | 828/852    | 10286/10324 | 3187/3283  |

~~构建卡方检验表：AS与isform的保守性有无关系~~

|        | 保守的isform | 不保守的isform |
| ------ | --------- | ---------- |
| 发生AS   |           |            |
| 没有发生AS |           |            |

```r
> table=matrix(c(11423,1445,16448,7232),nrow=2)
> chisq.test(table)
```

在多倍化后At、Dt亚组保守的isform数减少了，统计在两个二倍体多倍化后四倍体内亚组间的差异

1. 在二倍体中保守的转录本
   1. 只在At中存在
   2. 只在Dt中存在
   3. At、Dt也存在
2. 在二倍体中不保守的
   1. 只在At中存在
   2. 只在Dt中存在
   3. At、Dt都存在

在`category`目录中放了这几类

```bash
 python ~/scripte/Alternative/module/FEST3/fourHomologIsform.py  -c1 ../A2_vs_D5/conserve_isform.txt  -c2 ../A2_vs_At/conserve_isform.txt  -c3 ../At_vs_Dt/conserve_isform.txt  -c4 ../D5_vs_Dt/conserve_isform.txt  -u1 ../A2_vs_D5/A2_uniq.txt  -u2 ../A2_vs_D5/D5_uniq.txt -o1 1 -o2 2 -o3 3 -o4 4 -o5 5  -A2 A2_uniq -D5 D5_uniq
```

#### 比较保守的isform之间表达量是否存在差异

```bash
##脚本内对两个基因组总的表达水平进行了矫正
python ~/scripte/Alternative/module/FEST3/conserveIsformFPKM.py -c conserve_isform.txt  -fpkm1 ../../geneExpress/stringtie/A2/t_data.ctab  -fpkm2 ../../geneExpress/stringtie/TM1/t_data.ctab  -o hhhh
##获取差异倍数变化
awk -v fild1=5 -v fild2=6 -f ~/github/zpliuCode/script/log2.awk hhhh >hhhh2
```

![Nla5uR.png](https://s1.ax1x.com/2020/06/20/Nla5uR.png)

#### isform水平的差异表达

| 比较    | 总数    | 上调表达 | 下调表达 |
| ----- | ----- | ---- | ---- |
| A2 D5 | 11823 | 2341 | 2650 |
| At Dt | 6331  | 1114 | 1193 |
| A2 At | 12868 | 3273 | 2300 |
| D5 Dt | 12335 | 3182 | 2120 |

#### isform表达水平的剂量效应

#### 保守的isform，与不保守的isform平均表达水平差异

#### 特异性的isform

```bash
python ~/scripte/Alternative/module/FEST3/uniqisform.py -c conserve_isform.txt  -a isform_count.txt  -u1 A2_uniq.txt  -u2 D5_uniq.txt
```

### 参考

1. [wu-blast](https://www.advbiocomp.com/blast/obsolete/)&#x20;
2. [ab-blast](https://blast.advbiocomp.com/licensing/personal.html)


---

# Agent Instructions: Querying This Documentation

If you need additional information that is not directly available in this page, you can query the documentation dynamically by asking a question.

Perform an HTTP GET request on the current page URL with the `ask` query parameter:

```
GET https://zpliu.gitbook.io/booknote/ke-bian-jian-qie/di-san-ci-dui-as-jin-hang-tong-ji/02.md?ask=<question>
```

The question should be specific, self-contained, and written in natural language.
The response will contain a direct answer to the question and relevant excerpts and sources from the documentation.

Use this mechanism when the answer is not explicitly present in the current page, you need clarification or additional context, or you want to retrieve related documentation sections.