多倍化过程中变化2

保守转录本与不保守转录本的差异

特异性的转录本：A2与At相比CDS不保守的转录本

1.保守转录本与不保守转录本在亚基因组间的占比

多倍化过程中有54.1%~60.7%的转录本是特异性的转录本，并且在A基因组中有更高比例的特异性转录本存在，D5基因组中特异性转录本的比例最低

基因组	保守转录本	特异性转录本	总的转录本
A2	17915(39.3%)	27629(60.7%)	45544
At	13386(41.8%)	18612(58.2%)	31998
Dt	15357(45.9%)	18106(54.1%)	33463
D5	17151(44.0%)	21821(56%)	38,972

2.转录出特异性转录本的基因数，在多倍化过程中逐渐减少，多倍化后更少的基因转录出多个特异性的转录本

在A2中一共有8557个基因存在特异性的转录本，At中一共有7630个基因存在特异性的转录本，D5中一共有8103的基因存在特异性的转录本，Dt中一共有7420个基因存在特异性的转录本。

在二倍体中更多的基因转录出多个特异性的转录本

cut -f1 Dt_specialIsoformAnnotion.txt |sort |uniq -c|awk '{print $1}'|sort |uniq -c|awk '$2>=10{a+=$1}$2<10{print "Dt\t"$2"\t"$1}END{print "Dt\t>10\t"a}'

3.特异性转录本与保守转录本在CDS长度和FPKM上的比较

低表达的转录本：FPKM值<1，卡方测验

不同类型的转录本与转录本的表达水平是否存在差异

At中类型	低表达	高表达
特异性的转录本	10424(56.0%)	8188
保守的转录本本	5163(38.6%)	8223

Dt中类型	低表达	高表达
特异性的转录本	10519(58.0%)	7587
保守的转录本本	5854(38.1%)	9503

A2中类型	低表达	高表达
特异性的转录本	16706(60.4%)	10923
保守的转录本本	7572(42.3%)	10343

D5中类型	低表达	高表达
特异性的转录本	11935(54.7)	9886
保守的转录本本	6298(36.7%)	10853

##进行卡方测验
data=matrix(c(10424,5163,8188,8223),nrow=2)
chisq.test(data)
#或者fisher.test(data)

4.特异性转录本中不同AS类型的数目

在同一个直系基因组中将特异性的转录本与保守的转录本进行比较，发现差异的AS在其中的比例，

并且这种差异的AS被认为是亚组特异性的AS

基因组	保守的转录本	特异的转录本	特异转录本与保守转录本间存在AS 差异	其他原因导致的特异性转录本
A2	17915	27629	17566	10063
At	13386	18612	12831	5781
D5	17151	21821	15032	6789
Dt	15357	18106	12530	5576

这些AS事件的种类和数目的统计

##剪切事件数目的统计
cut -f3 A2_specialIsoformAnnotion.txt |sed 's/,/\n/g'|sort |uniq 
##存在AS事件差异数目，和特异性转录本数目之间的相关性
python ../../D5_vs_DT/conserve_vs_specialIsoform/special_AS.py A2_specialIsoformAnnotion.txt  11

同时基因产生特异性转录本的数目与AS的数目存在正相关，这只筛选存在AS差异的基于进行相关性分析

##相关性分析函数
cor.test()

基因组	IR	ES	A5	A3	total
A2	7231	1132	2788	4139	15290
At	5953	1133	3120	4044	14251
D5	6805	1452	3186	5046	16489
Dt	5979	1076	3003	3905	13963

5.A、D亚组间特异性转录本的丢失与获得

A2 vs At 11292

D5 vs Dt 11690

不同的基因组对于特异性转录本的比例是否存在相关性，在多倍化过程中A基因组显著性的更多的基因发生了特异性转录本的丢失；而在四倍体中同样的是At中更多基因获得了特异的转录本。

卡方测验：特异性转录本的丢失

A2 和D5进行比较谁丢的多的时候，A2和D5的总的基因数目要一样；

At和Dt进行比较谁获得的特异性转录本多的时候，At和Dt总的基因数目要一样多；

使用四元同源基因对。

基因组	存在特异性转录本的基因	不存在特异性转录本的基因	Total
A2	8557(75.77%)	2735	11292
D5	8103(69.3%)	3587	11690
At	7630(67.5%)	3662	11292
Dt	7420(63.5%)	4270	11690

基因组	存在特异性转录本的基因	不存在特异性转录本的基因	Total
A2	6927(64.44%)	3824	10751
D5	6519(60.63%)	4232	10751
At	6154(57.24%)	4597	10751
Dt	5984(55.66%)	4767	10751

awk -F "\t" '{split($2,a,",");split($5,b,",");if(length(a)==1&&length(b)==1){print $0}}' A2_special_At_allConserve.FPKM |less

特异性转录本丢失的例子

特异性转录本获得的例子

5.根据基因包含特异性转录本的情况对基因进行分类

1. high

2. mild

3. none

不同类别基因的表达水平的差异

6.基因主要表达的转录本是否发生改变

在多倍化过程中存在特异性转录本和保守转录本的基因中，基因所表达的productive 转录本是否发生了改变。

##提取既存在保守转录本又存在特异性转录本的基因对
awk -F "\t" '$4!=""{print $0}' A2_allConserve_At_special.txt |cut -f1-4  >Isoform_DomainFPKMChange/conserve_specialIsoform.txt 
awk -F "\t" '$2!=""{print $0}' A2_special_At_allConserve.txt |cut -f1,2,4,5 >>Isoform_DomainFPKMChange/conserve_specialIsoform.txt 
awk -F "\t"  '$5!=""&&$2!=""{print $0}' gene_conserve_allSpecial.txt |cut -f1,2,4,5 >>Isoform_DomainFPKMChange/conserve_specialIsoform.txt 
##这些基因的主要表达的转录本是否仍旧保守
python /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/evolution4/A2_vs_AT/conserve_vs_specialIsoform/special_vs_conserveFPKMDomain.py -Aisoform  ../A2_PacBio.txt   -Bisoform ../TM1_PacBio.txt  -input Isoform_DomainFPKMChange/conserve_specialIsoform.txt -o Isoform_DomainFPKMChange/conserve_specialIsoformDomain.txt

类型	主要表达的转录本发生改变	productive 转录本没有发生改变	Total
A基因组	1599	3338	4937
D基因组	1375	3604	4979

可以做一个胶图，那些有些高表达基因既存在保守转录本，又存在不保守转录本，但是表达量高的是不保守转录本

对这些主要表达转录本发生改变的基因进行功能富集分析

• A基因组中主要表达的转录本发生改变 http://kobas.cbi.pku.edu.cn/kobas3/retrieve/?taskid=06a7875bd08d45d9ab659b6ec67221a8
• D基因组中主要表达的转录本发生改变http://kobas.cbi.pku.edu.cn/kobas3/retrieve/?taskid=ff13d4ec8e3e4c538b841fba6c8830d7

处理GO的输出结果

##去除基因id信息
grep -E  -v "^[#-]" AgenomeGO.txt |sed '/^$/d'|cut -f1-7|awk -F "\t" '$7<=0.01{print $0}'|less
##找单个基因的注释信息
grep Ghir_A07G023100 /data/cotton/zhenpingliu/genome/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_PrimaryPeptide_vs_Arabidopsis_Annotation.txt

A基因组中基因富集在脂质代谢作用上；

D基因组中基因富集在激动蛋白途径，与细胞骨架相关

找差异表达基因做胶图

筛选发生剪接事件的地方有较多的read覆盖，同时基因的表达水平也比较高

##筛选剪接事件出比对的read count
python /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/16example/stat_ASevent_readCount.py  -a A2_AS.bed -b hisat2.bed -chromsome chromsome.bed -o AS_AS_count.txt
##

举例子

A2 vs At中的基因

evm.TU.Ga05G1374 vs Ghir_A05G013330 二倍体中高表达的转录本包含一个IR事件，而四倍体At中是正常的那个

evm.TU.Ga03G0628 vs Ghir_A02G005750 二倍体中存在一个高表达的IR转录本

evm.TU.Ga03G1862 vs Ghir_A03G015520 二倍体中存在一个高表达的IR转录本

evm.TU.Ga05G0904 vs Ghir_A05G008720 二倍体中存在一个高表达的A5转录本

evm.TU.Ga05G1931 vs Ghir_A05G018570 四倍体中存在一个RI高表达基因

PB.10935.3 、PB.11193.4

D5 vs Dt中的基因

Gorai.006G171900 vs Ghir_D09G015850 ,在四倍体中存在一个RI的高表达转录本

Gorai.003G002800 vs Ghir_D03G000260,在四倍体中存在一个RI高表达的转录本

Gorai.001G174600 vs Ghir_D07G015890,在四倍体中存在一个RI的高表达的转录本

7.同一个基因组内，不同类型基因的比较

多倍化过程中，稳定的转录后调控基因表达水平相比于不稳定的转录后调控基因表达水平更高；可能不稳定的转录后调控基因产生的转录本大多是无意义的

不存在保守转录本的这些基因，表达量最低，这里把两个基因组的数据放在同一个图内，并且对不同基因组间的表达量进行一个标准化

A2 vs At 时，保证基因的总的表达量是一致的

A2基因总表达量为：886732 read数为 77490057

At基因总表达量为：270235 read数为 55815433

Dt基因总表达量为：286558 read数为 58665229

D5基因总的表达量为：585450 read数为 73759547

At按照A2总表达量进行标准化，乘以3倍；Dt按照D5的表达量进行标准化乘以2

A、D两个亚基因组间的差异

大多数AS可能是noise，并且产生截断的蛋白质；并且物种间AS的不保守可能也是由于AS与物种特异性功能的形成有关；例如人类和老鼠中仅仅只有10%的AS event在同源基因间是保守的；AS的快速分化有助于物种的形成。

类型	数目
所有转录本都保守	1170
At中所有转录本保守、但是Dt中有特异性的转录本	855
Dt中有特异性的转录本、但是At所有转录本都保守	943
At、Dt除了保守转录本，两个都存在特异性转录本	1621
不存在保守的转录本	5062

根据转录本的保守情况对两个亚基因组同源基因进行分类：

• 存在保守转录本的同源基因

• 不存在保守转录本的同源基因

比较两类基因在二倍体和多倍体中的状态（比较7605组四元同源基因对，FPKM都>1,PacBio数目>1）：

 ##二倍体中存在保守转录本的基因
 evolution4/A2_D5_At_Dt/diploid_conserve_geneList

type	存在保守转录本	不存在保守转录本	total
二倍体	5157(67.8%)	2448(32.2%)	7605
四倍体	3018(39.7%)	4587（60.3%）	7605

多倍化后更多的亚基因组同源基因不存在保守的转录本；多倍化加深了A、D亚组间转录本的差异

• A2、D5中保守的基因在多倍化后，变得不保守

  • At、Dt的比例

• A2、D5中不保守的基因在多倍化后，变得保守。

  • At、Dt的比例

并行分化

二倍体中保守的基因，但在多倍化后At、Dt间发生了分化，分化的情况统计；

• 两个亚基因组与二倍体都不保守

• 其中一个亚基因组与二倍体存在保守的转录本，而另一个亚基因组与二倍体不存在保守的转录本

在5157个二倍体中保守的转录本中，其中有2670个在多倍化后变得不保守；在这些多倍化过程中变化的四元同源基因对中

python /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/evolution4/A2_D5_At_Dt/differentitation/differentitation.py -homolog differentiation_geneList.txt -A2D5 ../diploid_conserve.txt  -A2At ../A2_vs_At_conserve.txt  -D5Dt ../D5_vs_Dt_conserve.txt  -o test

At发生了分化，Dt与二倍体保守的转录本一样	1018(38.12%)	2670
类型	基因数目	total
Dt发生了分化，At与二倍体保守的转录本一样	961(36%)	2670
At、Dt发生分化，都与二倍体不保守	691(25.88%)	2670

通过卡方检验发现，At、Dt在多倍化过程中发生分化的基因数目上没有显著性差异，

比较这些基因的表达水平和转录本的表达水平；为了进行四组同源基因间基因表达水平的比较，将所有的同源基因按照A2的水平进行标准化

At中所有的基因表达水平乘以3.28，Dt中所有基因的表达水平乘以3.1，D5中所有基因的表达水平乘以1.51

A基因组发生分化的基因相比于D基因有着更低的表达水平；

D基因组发生分化的基因相比于A基因有着更低的表达水平

A、D基因组都发生分化的基因有着最低的表达水平。

多倍化后同源基因间转录本的分化表明它们剪切模式可能发生了改变

在At发生分化的基因中，有多少基因存在特异性AS，鉴定不同同源基因间的AS模式的差异。

类型	存在AS基因	特异性的AS的基因	total
A基因组发生分化	648	559	1018
D基因组发生分化	606	519	961

At发生分化的基因，特异性的AS事件

在At分化的基因，找到对应的特异性转录本，统计特异性转录本所包含的AS 事件，过滤掉那些与其他基因组保守的事件，最终得到的就是导致A基因组发生分化的AS 事件。亚组间特异的AS事件，找到它在其他基因组的位置，使用k-mer运算，在blast后得到相似的k-mer片段，在相似位置 处的k-mer则认为是对应的序列

对得到的k-mer序列与AS事件序列进行blastn，相似度大于90%，相似长度占AS事件长度达到90%，e-value<1e-5

##统计各个基因组中保守的AS 事件
awk '$2!=""{print $1 "\t"$2}' /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/evolution4/A2_D5_At_Dt/blast/At_Dt/At_Dt_conerved_AS >>At_conserve_AS.txt 
##发生AS的分化基因
cut -f3 ../A2_D5_Dt_conserve.txt |xargs -I {} grep {} ../../../../evolution2/TM1_AS.txt  >111
grep "PB" 111  >At_differentiation_AS.txt

Dt发生分化的基因，AS模式的改变

趋同进化

二倍体内转录本不保守，而在多倍化后At、Dt亚基因组存在保守的转录本；

• A2向D基因组方向趋同

• D5向A基因组方向趋同

##二倍体内转录本差异的基因
cat ../diploid_conserve_geneList ../filter_homologGene.txt |sort |uniq -u  >diploid_differentiation.txt 
##差异的基因在多倍化后变得保守
cat diploid_differentiation.txt  ../ployploid_conserve_geneList |sort |uniq -d >parallel_geneList
##A基因组发生趋同
python parallel.py  -A2 ../../A2_PacBio.txt  -D5 ../../D5_PacBio.txt  -TM1 ../../TM1_PacBio.txt  -ho parallel_geneList  -oA A2_At_Dt_geneList -oD D5_At_Dt_geneList
##D基因组发生趋同

二倍体中差异的基因，在多倍化后变成了保守的状态；

在2448个差异的二倍体基因中，共有531个基因在多倍化后发生了趋同进化，其中有174个基因发生了D5向A2方向的趋同，有266个基因发生了A2向D5方向的趋同，仍旧有99个基因趋同的方向是不确定的；

多倍化过程中，更多的同源基因向D基因组方向同化

类型	基因数目	total
向A基因组同化	174（32.8%）	531
向D基因组同化	266（50.1%）	531
无方向的趋同	99（17.1%）	531

At、DT同源基因间AS的分化

两个亚基因组同源基因间的AS是否发生了分化，分析每个基因组同源基因所包含的AS数目，以及AS保守的比例。

以下为21066个四元同源基因，所包含剪切事件的数目和基因数；

这里的AS 事件只使用发生在PacBio转录本上的

同源基因	剪切事件数目	基因数
A2	18814	6417
D5	18991	6805
At	13981	5122
Dt	14083	5115

使用k-mer获取对应保守的AS事件

 python conserve_AS_kmer.py  -A2 A2_AS.txt  -At At_AS.txt  -Dt Dt_AS.txt  -D5 D5_AS.txt  -kmer A2_kmer.txt  -o 111
##获取各个基因组间保守的剪切事件
awk '$2~/Ghir_A/{print $1"\t"$2}$3~/Ghir_A/{print $1"\t"$3}$4~/Ghir_A/{print $1"\t"$4}' A2_conserve_kmer_AS.txt

awk '$2~/evm/{print $2"\t"$1}$3~/evm/{print $3"\t"$1}$4~/evm/{print $4"\t"$1}' At_conserve_kmer_AS.txt

## 存在保守转录本的基因中，保守AS模式的基因
cut -f3 ployploid_conserve_geneList|xargs  -I {} grep {} At_Dt_AS_diversity.txt |awk '$3=="noDiver"{print $0}'|wc -l

同源基因间存在AS的差异

AS差异：存在不保守的AS事件

对这7605个同源基因AS的保守性分析，发现有5032（3166+1866）对基因在AS上存在差异；并且不存在保守转录本的同源基因对中AS差异的比例，显著的高于存在保守转录本的同源基因对

类型	AS差异	AS不存在差异	Total
存在保守转录本	1866	1152	3018
不存在保守转录本	3166	1421	4587

在不存在保守转录本的基因中，AS的差异显著的高于存在保守转录本；说明亚基因组间AS的差异是导致两个基因组产生不同转录本的原因之一。这里举个例子，由于AS的分化导致At、Dt转录本的差异

有多少同源基因的分化是由于AS导致的转录本的分化

A2、D5存在保守转录本、At与Dt间不存在保守转录本；并且其中一个基因组是由于特异性的AS导致

接下来主要讨论A、D亚组间productive转录本的差异

productive转录本在多倍化过程中的变化

在CDS的长度是否发生改变：

• A2和D5中productive 转录本一样；At和Dt中存在一个一样的productive 转录本（单个受AS调控）

  At分化的基因中的productive 转录本，与参考转录本相比是否存在AS

• A2和D5中productive 转录本一样；At和Dt中都不存在一样的转录本（两个基因同时受AS调控）

• A2和D5中productive 转录本不一样，但是At与Dt分别与对应的祖先一样（多倍化前就发生了分化）

• A2和D5中productive 转录本不一样，多倍化后只保留了A基因组或者只保留了D基因组的productive 转录本（单个基因受到AS调控）

• A2和D5中productive转录本不一样，At和Dt也都找不到与祖先一样的（同时受到AS调控）

在7605对At、Dt分化的基因中

#统计分化基因的produc转录本
 python geneDivesityByAS.py   -TM1 ../../evolution4/TM1_PacBio.txt  -A2 ../../evolution4/A2_PacBio.txt  -D5 ../../evolution4/D5_PacBio.txt -homo ployploid_nonConserve_geneList  -o 11

即使存在保守的转录本，而productive 转录本仍旧可能发生改变

类别	数目
At、A2、D5	754
Dt、A2、D5	822
只有A2、D5相同的production 转录本	674
A2、D5不一样，继承给了At、Dt	257
A2、D5不一样，At与A2一样，Dt与D5仍旧不一样	662
A2、D5不一样，Dt与D5一样，At与A2仍旧不一样	947
A2、D5不一样，At与Dt也都没找到与祖先一样的	2132
A2、D5、At、Dt都一样	1357
total	7605

在7605个四元同源基因对中，将两个亚基因组与二倍体祖先基因组相比。

• Dt亚组保守而At亚组转录本发生变化数目（822+947=1769）

• At亚组保守而Dt亚组转录本发生变化数目（754+662=1416）

• At和Dt亚组转录本与祖先相比都发生变化（674+2132=2806）

• At和Dt都没有变（1357+257=1614）

多倍化加剧了At、Dt亚组转录本的分化，将四倍体基因与二倍体祖先相比发现2806（36.89%）的亚基因组同源基因productive isoform同时发生了分化，有1769（23.26%）At同源基因的productive isoform发生了分化；而有1416（18.62%）的Dt同源基因productive isoform发生了分化；存在1614（21.22%）的同源基因在多倍化过程中保持一致的productive isoform。

在At与Dt之间仅有2088（27.46%%）的基因存在保守的productive isoform；但是（1357+947+662+257+822+754）63.11%的基因在多倍化过程中没有改变productive isoform。功能冗余的基因组，在同一个基因组后由于其中一个基因组受到的选择压力变小了，从而发生了分化。

总共（1769+1416）个同源基因发生了production的分化；更多的At基因在productive转录本上发生了分化;

比较At发生分化的productive 转录本FPKM与Dt中没有发生分化的productive转录本的FPKM的差异；

###获取转录本表达量的log2差异倍数
awk '$7<1&&$8<1{print 1}$7<1&&$8>1{print 1/$8}$7>1&&$8<1{print $7/1}$7>1&&$8>1{print $7/$8}' 11

在At发生分化的基因中，At的转录本表达水平显著性的低于保守的Dt中的转录本；

在Dt发生分化的基因中，Dt转录本的表达水平显著性的低于保守的At；在At亚基因组中productive isoform发生分化的基因中，这些分化产生的isoform相比于保守的productive显著的发生了截断。

大部分发生分化的productive isoform可能由于AS引入提前终止的密码子，导致表达水平的下调和CDS序列的截断。

韦恩图，统计At与Dt间productive isoform的保守情况。

有多少gene Productive isoform分化的可能是由于AS产生的

• 将Production isoform发生改变的住转录本与参考转录本进行AS的比较。

python /public/home/zpliu/work/Alternative/result/Gh_result/CO31_32_result/evolution5/At_Dt/ASregulation/At_differentation/ASregulateproductiveIsoform.py  -ref ../../../../evolution4/TM1_reference_isoform.txt  -PacBio ../../../../evolution4/TM1_PacBio.txt  -AS ../../../../evolution2/TM1_AS.txt  -homo A2_At_productiveIsoform.txt  -o 11

At中production isoform发生改变；而Dt中production 转录本是保守的；鉴定到了202个基因322个AS事件

Dt中production isoform发生改变；而At中production转录本是保守的；鉴定到了186个基因283个AS事件

总共鉴定到了388个基因对应605个剪切事件可能，导致production isoform 发生改变；

这605个剪切事件中有多少比例是亚基因组特异性的剪切事件。

举个特例、跑个胶图

At production isoform 发生改变的基因中；

Ghir_A13G024480存在一个IR事件，导致产生一个新的productive转录本；

Ghir_A11G008990存在一个特异的IR事件，而Dt中不存在这个IR事件

A2、D5、At、Dt四个AS模式是否保守：

• A2、D5、At（At基因组中保守的剪接模式，亚基因组间AS的分化）

• A2、D5、Dt（Dt基因组中保守的剪接模式，亚基因组间AS的分化）

• A2、At、Dt（AS的趋同）

• D5、At、Dt（AS的趋同）

• A2、D5、At、Dt（四个基因组都保守的模式

对于不保守的事件，使用k-mer找到对应的坐标