第7个结果

• 各个剪切事件并不是显著性的差异

• ES数据核实一下

• 同源基因之间的poly A的差异

• 组成型intron与剪切intron之间长度差异

• 同源基因在各个棉种的比较中AS的差异

• 同源基因间表达的差异与保守度

• 各个基因组中TE的插入

处理A2的原始数据

参考

https://www.omicsclass.com/article/344

AS数据的处理，由于D5基因组相比A2基因组，注释的转录本更多，这里看一下只考虑PacBio测序得到的与参考基因组进行比较的结果

for i in A2 D5 TM-1; do for j in IntronR AltA AltD AltP ExonS Other; do printf ${j}"\t"; grep ${j} ../${i}/end_third|grep PB|wc -l; done; done

和对应的基因数目

for i in A2 D5 TM-1; do for j in IntronR AltA AltD AltP ExonS Other; do printf ${j}"\t"; grep ${j} ../${i}/end_third|grep PB|cut -f2|sort|uniq|wc -l; done; done

基因组	IR	AltA	AltD	AltP	ES	Other
A2	28078	5199	4135	1823	2276	1435
D5	22226	3257	2434	1131	1719	1039
TM-1	40354	7823	6854	4665	3071	2539
At	19869	3868	3341	2301	1531	1197
Dt	20436	3946	3505	2363	1540	1342

看看公司跑的数据，最后就使用这个数据了

鉴定可变剪切的原理:

1. 将PacBio测得的isoform比对到参考基因组，获得isform 注释信息

2. 将PacBio isoform信息与参考基因组组中原有的注释信息进行合并

3. 使用脚本，根据合并后的注释信息中的基因结构信息，提取对应剪切事件和发生的位置

由于A2基因组中原有的注释信息里没有不同转录本的注释信息，并且PacBio测序只测了leaf这一个组织，所以在与A2中原有的注释信息合并后，也不能完全的代表所有的转录本的注释信息

而D5和TM-1基因组中有不同转录本的注释信息，因此在根据基因结构去鉴定AS的时候，有有一些差异。

为了减少这个带来的差异，我在各个棉种中选择那种只被PacBio isoform支持的AS

没有统计scaffold上的AS事件

• AS事件数目

基因组	IR	AltA	AltD	AltP	ES	Other
A2	36574	5847	4472	4394	2352	4269
D5	30379	5350	3845	4185	2229	6302
TM-1	51453	9749	8252	7358	3723	7915
At	25386	4878	4011	3621	1850	3882
Dt	26067	4871	4241	3737	1873	4033

• gene 数目

基因组	IR	AltA	AltD	AltP	ES	Other
A2	9139	3147	2614	2200	1499	1881
D5	7466	2900	2295	2214	1443	2285
TM1	13574	5414	4740	3780	2399	3308
At	6657	2709	2345	1872	1210	1672
Dt	6917	2705	2395	1908	1189	1636

提取FESTs序列

python ~/scripte/Alternative/module/extractFEST2.py -p ~/work/Alternative/result/Gh_result/CO31_32_result/06_Alignment/all.collapsed.gtf  -r ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf -AS ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt  -g ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta  -IR 1 -ES 2 -AltA 3 -AltD 4
## D5在分析的时候，取的是CDS坐标

##拆分At与Dt

blast分析保守性

• 相似度80% ，覆盖度80%以上 e-value 1e-5

统计每种AS对应的比例

# evm.TU.Ga01G0013        0       Ghir_A01G000130 0
出现这种情况是因为，这个基因不存在IR事件，但是存在AltA获取AltD等其他事件
## 统计同时存在AS的基因数，和AS数目
for i  in `ls .`; do  printf ${i}"\t"; awk '$2!=0&&$4!=0{a+=1;b+=$2+$3}END{print a"\t"b}' ${i}; done
## 统计保守的基因数，与AS数

保守IR事件

存在保守事件基因的比例：

再计算每个基因，AS事件的保守程度：

保守基因的比例:

1. 2425/7142 33.95%

2. 2274/6594 34.48%

3. 31.91%

4. 62.59%

保守事件的比例:

1. 4789/36144 *2 26.50%

2. 26.95%

3. 24.14%

4. 50.14%

比较	基因对数	事件数	保守基因对数	保守事件数没乘2
A2 vs At	7142	36144	2425	4789
D5 vs Dt	6594	32642	2274	4398
A2 vs D5	7218	35364	2303	4269
At vs Dt	2911	13551	1822	3397

保守ES事件

保守基因的比例:

1. 31.10%

2. 32.11%

3. 23.28%

4. 56.40%

保守事件的比例:

1. 33.85%

2. 36.47%

3. 26.42%

4. 65.40%

比较	基因对数	事件数	保守基因对数	保守事件数没乘2
A2 vs At	598	1152	186	195
D5 vs Dt	570	1146	183	209
A2 vs D5	524	984	122	130
At vs Dt	250	474	141	155

保守AltA事件

保守基因的比例：

1. 55.22%

2. 48.02%

3. 41.05%

4. 45.61%

保守事件的比例:

1. 48.47%

2. 48.25%

3. 41.08%

4. 45.98%

比较	基因对数	事件数	保守基因对数	保守事件数没乘2
A2 vs At	891	2030	433	492
D5 vs Dt	835	1919	401	463
A2 vs D5	765	1704	314	350
At vs Dt	820	1792	374	412

保守AltD事件

保守基因的比例：

1. 49.92%

2. 46.88%

3. 40.26%

4. 46.22%

保守事件的比例:

1. 52.32%

2. 51.19%

3. 44.55%

4. 49.36%

比较	基因对数	事件数	保守基因对数	保守事件数没乘2
A2 vs At	659	1380	329	361
D5 vs Dt	640	1344	300	344
A2 vs D5	544	1055	219	235
At vs Dt	662	1414	306	349

比较不同的AS类型，保守程度是否存在差异

• 计算保守基因对，的AS事件保守程度

• 不同剪切事件在多倍化的过程中，保守的程度不一样，其中ES最保守，而IR保守性最低

不同AS事件，保守性不一样

• 同一类剪切事件，在不同基因组的比较中的保守程度差异

A2 vs At 叫保守，D5 vs Dt叫保守；

A2 vs D5叫并行分化，At vs Dt叫趋同进化

不同剪切事件的长度差异

• 组成型的内含子

• 保守的剪切事件对应的intronR

• 不保守的剪切事件对应的intronR

#所有的型内含子
extract_splice_sites.py ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf  >1111
awk '{print $1,$2+2,$3,$4}' OFS="\t"  1111 >all_intron.bed
 awk '{print $1,$2+2,$3,$4}' OFS="\t" ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model_splce.txt >> all_intron.bed
sort -k1,1 -k2,3n all_intron.bed |uniq >1111
mv 1111 all_intron.bed

 ~/software/bedtools2-2.29.0/bin/intersectBed  -loj -a ../all_intron.bed  -b ./As_intronR.bed  >222
# 发生剪切的内含子
grep ">" ../../../TM1/1|sed 's/>//g'|awk -F "-" '{print $2"\t"$3"\t"$4"\tAlternative"}' |uniq  >As_intronR.bed

#组成型的intron
 awk '$5=="."{print $0}$5!="."&&($7-$6)/($3-$2)<0.4{print $0}' 222|cut -f1,2,3|sort |uniq |awk '{print $0"\tConstitutive"}'>constitutive_intron.bed
#存在剪切的内含子中保守的
cut -f3 ../../../blast/IR/A2_At_conserve_end|awk -F "-" '{print $2"\t"$3"\t"$4"\tconserve"}' >conserveAS.bed
#存在剪切的内含子中不保守的
 cat conserveAS.bed  As_intronR.bed |sort |uniq -u >noconserveAS.bed

• 组成型的外显子

• 发生ES的外显子

• 保守的ES

##所有exon坐标
 cat ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf |awk '$3~/exon/{print $1,$4,$5,$7}' OFS="\t" >all_exon.bed
 sort -k1,1 -k2,3 all_exon.bed |uniq >111
 mv 111 all_exon.bed
 ## 发生剪切的exon的坐标
 cat ../../../TM1/ES_length >./As.bed

 ##组成性exon
 cat all_exon.bed  As.bed |cut -f1,2,3|sort |uniq -u |awk '{print $0"\tconstitutive"}' >constitutive.bed

 ## 保守的ES的坐标
cut -f3 ../../../blast/ES/A2_At_conserve_end|xargs  -I {} grep {} As.bed |awk '{print $0"\tconserve"}' >conserve.bedd
 ## 不保守的ES坐标
cut -f1-4 conserve.bed|cat -  As.bed |sort |uniq -u |awk '{print $0"\tnocnoserve"}' >noconserve.bed

对多倍化过程中的剪切事件进行分类

• parallel conserve

  • A2与D5中保守，在At与Dt中同样保守

  • A2与D5中不保守，但是在At与Dt中保守；A2 At或者D5 Dt中存在保守

• bias conserve

  • A2与D5中保守，但是在At与Dt中不保守；A2 At或者D5 Dt中存在保守

  • A2与D5中不保守，在At与Dt中同样不保守；但是在A2 At或者D5 Dt中是保守的

## 目录 ParallelEvolution

四倍体内同源基因更多的发生分化

分析A2与D5中保守的剪切事件，在At、Dt中的保守情况 1. A2 D5中保守，在At与Dt中同样保守 2. A2 D5中保守，在At与Dt中不保守

A2和D5基因组分别是一个独立的个体，存在一些保守的AS，维持正常的生命活动；

剪切机制的存在往往是为了让植物能够快速响应外界环境的变化，在正常的环境中这种机制的存在往往是能量的浪费，因此在四倍体中，发现大多数原来在二倍体内保守的剪切事件，在四倍体内发生了分化，只在At中存在，或者只在Dt中存在，这两种比例各占一半

A2和D5作为一个独立的个体，因此具有一些保守的剪切事件；而剪切事件往往是一种能量的浪费；在四倍体内只需要某一个基因组存在对应的AS即可

类型	二倍体内保守	四倍体内保守	四倍体内不保守
IR基因	2303	536	1154
IR事件数	4269	838	1697
ES基因	122	21	46
ES事件	130	22	47
AltA基因	314	48	121
AltA事件	350	51	127
AltD基因	219	37	77
AltD事件	235	40	80

绘制图形

在画一个维恩图，A中保守的，D中保守的，At、Dt中同样保守的比例

##A基因组中保守的事件数
cat 1 2 3 4 |grep "Ghir_A"|grep "evm"|awk '{print $1"^"$2"^"$3"^"$4}'>Agenome_conserve.txt
cat 1 2 3 4 |grep "Ghir_D"|grep "Gor"|awk '{print $1"^"$2"^"$3"^"$4}' >Dgenome_conserve.txt

四倍体内存在分化

事件数目统计

看剪切事件在两个亚基因组间是否存在偏向性

只在A基因组中保守的剪切事件，占A2和At中保守事件的比例

只在D基因组中保守的剪切事件

事件类型	只在A保守的占A中所有保守比例	只在Dt中保守的占D中所有保守的比例
IR	3623/4461	3189/4027
ES	164/176	154/186
AltA	394/445	358/409
AltD	287/327	270/310

剪切事件是否存在偏向性

分析每个基因组中特异剪切事件对应的GO

• 在A基因组中特异的AS对应的基因富集的GO

• 在D基因组中特异的AS对应的基因富集的GO

• A、D中保守的AS对应的基因富集的GO

#在A基因组特异的
cut -f3 1|awk '{print "^^"$0"^"}'|cat - Agenome_conserve.txt | awk -F "^" '{print $3}'|sort|uniq -u|cut -f1 -d "-"|sort|uniq |xargs  -I {} grep {} ~/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_Noscagenes_GO_V3_At.annot >A.GO
#在D基因组中特异的
cut -f4 1 |awk '{print "^^^"$0}'|cat - Dgenome_conserve.txt |awk -F "^" '{print $4}'|sort|uniq -u|cut -f1 -d "-"|sort|uniq |xargs  -I {} grep {} ~/genome_data/Ghirsutum_genome_HAU_v1.1/Gh_Noscagenes_GO_V3_Dt.annot >D.GO

分析可变剪切对基因表达式水平的影响

• At与Dt中存在保守事件的基因

• At与Dt中仅仅只有一个亚基因组中存在保守事件的基因

存在保守事件的同源基因对间的差异表达情况

不存在保守事件的同源基因对的差异表达情况

Divergence of duplicated genes by repeated partitioning of splice forms and subcellular localization

## 获取同源基因对
cat 1 2 |awk '{print $3"\n"$4}' |awk  -F "-" '{print $1}'|sort |uniq
cat 3 4 |awk '{print $3"\n"$4}' |sed '/*/d'|awk  -F "-" '{print $1}'|sort |uniq

## 获得同源基因对间的表达量
python homolog_expresion.py  -homolog ~/work/Alternative/result/homologo/homologGene/At_vs_Dt_collinerity.txt  -gene IR/At_Dt_bias_gene_ID  -FPKM ~/work/RNA-seq/hisat2_out/leaf/leaf_BAM/Gh/gene_fpkm.txt  -o IR/bias_fpkm

剪切事件在多倍化过程中的变化

• 在二倍体中存在的事件，并且在四倍体中At、Dt保守；多倍化之前就存在的事件

• 在二倍体中不存在的事件，但是在四倍体中保守；多倍化后新产生的事件

# 上个分类中1、2 文件

# 新产生的做个排除法

参考

方差分析 https://zhuanlan.zhihu.com/p/57756620

卡方检验 https://zhuanlan.zhihu.com/p/42803826