04

保守的AS

  • 使用wu-Blast进行对FEST进行比对

  • 首先提取转录本的注释信息

  • 根据AS发生的坐标提取FEST进行比对

  • 筛选两端都比对上的AS

python3 ~/scripte/Alternative/module/FEST3/conserveAS/conserveAS.py -p1 ~/work/Alternative/result/Ga_result/CO11_12_result/07_annotation/merge.gtf -p2 ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf  -ho ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt  -r1 ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gtf -r2 ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf -AS1 ~/work/Alternative/result/Ga_result/CO11_12_result/11_AS/end_splice.txt3  -AS2 ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt3  -g1 ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta -g2 ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -IR IR.txt -AltA AltA.txt -AltD AltD.txt -ES ES.txt
## 根据比对的长度，选择比对更长的
for i in `ls .`; do cat ${i} |awk '$1~/^Ghir_A/&&$2~/^Ghir_D/{print $0}$1~/^Ghir_D/&&$2~/^Ghir_A/{print $2,$1,$3}' OFS="\t"|sort |uniq |awk '{print $2,$3,$1}' OFS="\t"|sort  -k3 -k2,2nr|uniq -f2|awk '{print $3"\t"$2"\t"$1}'|sort -k3 -k2,2nr|uniq -f2  >${i}_end; done

## 统计同源基因间发生AS的数目
 for i in ExonS IntronR AltA AltD; do python ~/scripte/Alternative/module/homologASeventCount.py -homolog ~/work/Alternative/result/homologo/homologGene/D5_vs_Dt_collinearity.txt -AS1 ~/work/Alternative/result/Gr_result/CO41_42_result/11_AS/end_splice.txt3   -AS2 ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt3 -T ${i} -o ${i}; tmp=`mktemp`done; sort ${i}|uniq >${tmp}; mv ${tmp} ${i}; done

这个问题之后再看看

###这种情况得调整一下，两个IR事件，坐标实际上只相差了一个碱基；感觉可以合并后整成一个
-    -    Gorai.001G041400-Chr01-3801948-3802062    Ghir_D07G004050-Ghir_D07-4143773-4143886
-    -    Gorai.001G041400-Chr01-3801948-3802062    Ghir_D07G004050-Ghir_D07-4143773-4143887

保守的IR

比较

基因数

保守事件数

IR数1

IR数2

A2 At

2098

4786

12553

10749

D5 Dt

1611

4179

9642

8865

A2 D5

1936

4213

11552

10929

At Dt

1347

3171

7687

7491

保守ES

比较

基因数

保守事件数

ES数1

ES数2

A2 At

129

137

217

214

D5 Dt

111

122

211

224

A2 D5

98

107

164

161

At Dt

87

93

146

145

保守AltA

比较

基因数

保守事件数

AltA数1

AltA数2

A2 At

285

324

679

678

D5 Dt

222

253

564

524

A2 D5

211

236

488

513

At Dt

235

260

547

531

保守AltD

比较

基因数

保守事件数

AltD数1

AltD数2

A2 At

229

251

494

487

D5 Dt

162

183

403

386

A2 D5

140

155

264

305

At Dt

187

206

432

444

计算IR和ES事件的PIR值

  • 比较全基因组中所有intron的PIR值

  • 发生IR的IR的PIR值

保守事件对应的长度

TE insertion in retained introns may be a significant phenomenon in higher plants.

## 提取保守的IR事件对应的坐标
cat ../../conserveAS/A2_At/IR.txt_end ../../conserveAS/A2_D5/IR.txt_end|cut -f1 |sort |uniq |awk -F "-" '{print $2,$1,$3,$4}' OFS="\t" >conserve_IR
## 提取不保守的IR坐标
awk '$3~/IntronR/{print $1,$2,$4,$5}' OFS="\t" ../../../PSI/IR/A2_PSI.txt |cat - conserve_IR |sort -k1,2 -k2,3n|sort |uniq -u >noconserve_IR.txt
##提取Constitutive intron
+ 长度小于8000
+ PIR 小于0.1
awk '$3-$2+1<8000&&$6<0.1{print $0}' ../../../PSI/allIntron/A2_PSI.txt |awk '{a+=$3-$2+1}END{print a/NR}'

提取保守AS的PRI值

过滤掉所以那些基因中所有内含子区域都是为[0,0]的，这类可能是基因没有表达导致的

##没有表达的基因，从文件内过滤掉
awk '{a[$4][1]+=$5;a[$4][2]+=$6}END{for(i in a){if(a[i][1]!=0||a[i][2]!=0){print i}}}' TM1_PIR.txt |xargs  -I {} grep {} TM1_PIR.txt  >TM1_PIR2.txt

分析IR长度与PIR值的关系

  • 越短的IR，它的PIR值越高

模拟临界PIR值

从所有的intron中抽取1000个intron，用PIR阈值进行判断，得到200个为IR，而其中200个判断的IR有多少比例的是真正的IR，每个PIR处进行1000次迭代

随着PIR阈值的不断增加，真实的IR比例占潜在的比例反而越来越小；

  • 三代中测到了IR但是

PRI值增加，表明这个conserve AS多倍化后存在功能,conserve_isform_conserve_IR

cut -f1  ~/work/Alternative/result/homologo/FEST3/AS/conserveAS/A2_At/IR.txt_end|awk -F "-" '{print $2,$1,"IntronR",$3,$4}' OFS="\t"|xargs  -I {} grep {} ../IR/A2_PSI.txt >A2_conserve_IR.txt

AS保守性分析

  • A2与D5不保守，在四倍体中仍旧维持At与Dt不保守的状态，cDAS （conserve diversity AS）

  • A2与D5保守，在四倍体中At与Dt也保守ccAS（complete conserve AS）

  • A2、D5 At中存在AS，但At、Dt中只有一个存在AS BCAS（bias lost AS）

  • At与Dt同时丢失，At与Dt同时获得 contemporary gain AS | contemporary lost AS

  • Other

分类

A基因组

D基因组

CCAS都保守

1013

-

BLAS 偏向性的丢失

1358

1321

CDAS 保持AS的差异

571

908

CGAS 同时获得AS

1123+445+310

-

CLAS 同时丢失AS

556

-

Other

1357+946

-

找到所有同源基因在AS处对应的坐标

## 这种结果需要过滤一下，理论上是保守的IR，应该是取FEST的时候用300bp的差异导致的
Ghir_A01G004300-Ghir_A01-5344342-5345370    Ghir_D01G004400-Ghir_D01-4984060-4985096    y    n
Ghir_A01G004300-Ghir_A01-5344342-5345370    Ghir_D01G004400-Ghir_D01-4984060-4985095    n    y
## uniq掉
awk '{print $0"\t"$2}' A2_D5/A2_vs_D5_ASIR.txt |sort -k5,5 -k3,4r|uniq -f4|awk '{print $1,$2,$3,$4,$1}' OFS="\t" |sort -k5,5 -k3,4r|uniq -f4 >A2_D5/11111
##提取四组AS的情况
python ~/scripte/Alternative/module/homolog/identifyFourAS.py  -A2D5 ./A2_D5/11111  -AtDt ./At_Dt/11111 -A2At ./A2_At/11111  -D5Dt ./D5_Dt/11111  -o 2222
## 对应是否发生IR这个还得再用AS的数据看看
python judgeIR.py  -i A2_D5_At_Dt_AS.txt  -all all_IR.txt  -o 11

保守的AS和保守的isform

## 在鉴定的保守isoform中，同时存在多少保守的IR事件
python ~/scripte/Alternative/module/FEST3/conserveAS/isoform_ASconserve.py -As1 ~/work/Alternative/result/Ga_result/CO11_12_result/11_AS/end_splice.txt3  -As2 ~/work/Alternative/result/Gh_result/CO31_32_result/11_AS/end_splice.txt3  -isform ~/work/Alternative/result/homologo/FEST3/isforms/A2_vs_At/conserve_isform.txt -AS ../conserveAS/A2_At/IR.txt_end  -t IntronR -o 111

比较

保守isoform

保守的IR

A2 At

12868

481

D5 Dt

12335

327

A2 D5

11823

263

At Dt

6331

191

GO功能富集分析

将这些存在保守IR事件的四倍体基因提出来，进行GO富集

## At和Dt
## A2和At
## D5和Dt

特异性的AS在另外一个基因组的坐标

  • 提取一个基因组的所有intron坐标

##提取一个基因组的所有intron坐标
python3 ~/scripte/extractIntronbed.py  -fasta ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.0.fa -p ~/work/Alternative/result/Gr_result/CO41_42_result/07_annotation/merge.gtf  -r ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gtf -o 111 -f 2222
##去重
sort -k1,1 -k2,3n 111 |uniq|sed '/scaffold/d' >all_intron.bed

  • 提取特异性的AS坐标

    • 这里FEST序列就提取上下各300bp

mkdir A2_uniq
mkdir At_uniq
## 特异性AS坐标
awk '$3~/IntronR/{print $2"-"$1"-"$4"-"$5}' ~/work/Alternative/result/Ga_result/CO11_12_result/11_AS/end_splice.txt3|cat - ./IR.txt_end |cut -f1|sed 's/-[+-]$//g'|sort |uniq -u

##进行wu blast
python3 ~/scripte/Alternative/module/FEST3/conserveAS/noconserveAS.py -ho ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt  -s IR.txt  -all ../../../../PSI/allIntron/TM1_intron.txt -g1 ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta -g2 ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_genome_HAU_v1.0.fasta -o 111
##提取最保守的
cat 111|awk '$1~/^evm/&&$2~/^Ghir_A/{print $0}$1~/^Ghir_A/&&$2~/^evm/{print $2,$1,$3}' OFS="\t"|sort |uniq |awk '{print $2,$3,$1}' OFS="\t"|sort  -k3 -k2,2nr|uniq -f2|awk '{print $3"\t"$2"\t"$1}'|sort -k3 -k2,2nr|uniq -f2 >A2_IR_uniq_At.txt
##还得再过滤一下，因为all intron里可能不包含AS的那个intron
cut -f1 A2_uniq_D5_IR.txt |xargs  -I {} grep {} ../IR.txt_end |sed 's/-[+-]//g' |cat - A2_uniq_D5_IR.txt |sort -k1,1 |awk '{print $2"\t"$3"\t"$1}' |uniq  -f2 -u |awk '{print $3"\t"$1"\t"$2}' >111
cut -f3 D5_uniq_A2_IR.txt |xargs  -I {} grep {} ../IR.txt_end |cat - D5_uniq_A2_IR.txt |sort -k3,3 |uniq -f2 -u >111
mv 111 D5_uniq_A2_IR.txt
mv 111 A2_uniq_D5_IR.txt
##还有At uniq Dtuniq中的结果取重复就是保守的IR
cat D5_uniq_A2_IR.txt  ../A2_uniq/A2_uniq_D5_IR.txt |sort |uniq -d >>../IR.txt_end

不保守的AS的坐标统计

不保守的AS的原因

  1. 序列水平上的差异

  2. 表观上的差异

    sequence variation on DNA methylation analysis, we used the conserved sequences to examine DNA methylation changes between diploid and allotetraploid cottons.

甲基化差异与PIR差异

  • 看一下怎么鉴定甲基化差异基因的方法

519 differentially methylated genes identified between wild and cultivated cottons

100,246 CG, 109,424 CHG, and 252,042 CHH DMRs between allotetraploid and diploid species

在A2与D5中存在的DMRs ，在四倍体中仍旧维持的有20% cDMRs

A、D亚基因组在多倍化后发生反方向的变化，A到At减少，D到Dt增加；hDMRs；这有可能是染色体交叉互换的结果

Previous鉴定DRMs区域Next重新下载原始数据进行比对

Last updated