第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保守性最低

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

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

PreviousstringtieNext原始数据处理

Last updated