保守AS的鉴定
参考： Identification and analysis of alternative splicing events in Phaseolus vulgaris and Glycine max
Genome-wide identification of evolutionarily conserved alternative splicing events in flowering plan
提取两端exon序列
分别对应5种保守的ASs
其中AltD
2^,1^; 1^,2^;发生AltD的转录本为2开头的那个,就提取开头为1的转录本
AltA
2-,1-; 1-,2- AltA发生在转录本为1开头的那个，提取转录本为2的
IntronR
1^2-,0 ; 0,1^2-提取转录本为1^2-的那个转录本
Exon
1-2^3-4^,0 提取转录本为0的，并且ExonS指定的是包含那个exon的intron范围
##内部经过筛选除了，选择两端长度大于20bp的
python ~/scripte/Alternative/module/extractFEST.py
-p PAcBio.gtf
-r reference.gtf
-AS end_third
-g genome.fasta
-IR 1 -ES 2 -AltA 3 -AltD 4
##统计sequence Tag范围
grep -v ">" 1 |awk '{print length($0)}'|sort -n|less进行Blast
同一种类型的剪切、不同亚基因组之间建库
## A2 vs At 的intronR库
makeblastdb -dbtype nucl -in A2_At.fa -out A2_At
makeblastdb -dbtype nucl -in D5_Dt.fa -out D5_Dt
#区分At、Dt
grep "Ghir_A" -A1 4 |sed '/^--/d' >At4
##比对
blastn -query A2_At.fa -db A2_At -evalue 1e-5 -num_threads 2 -outfmt '6 qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore' -out A2_At.blast
blastn -query D5_Dt.fa -db D5_Dt -evalue 1e-5 -num_threads 2 -outfmt '6 qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore' -out D5_Dt.blast
sed -i 's/\.v2\.1//g' D5_Dt.blast
## 得到保守的AS
python ~/scripte/Alternative/module/homologBlast.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt -Blast A2_At.blast -o 222
python ~/scripte/Alternative/module/homologBlast.py -homolog ~/work/Alternative/result/homologo/homologGene/D5_vs_Dt_collinearity.txt -Blast D5_Dt.blast -o 222
python ~/scripte/Alternative/module/homologBlast.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_D5_collinearity.txt -Blast A2_D5.blast -o 222
python ~/scripte/Alternative/module/homologBlast.py -homolog ~/work/Alternative/result/homologo/homologGene/At_vs_Dt_collinerity.txt -Blast At_Dt.blast -o 222
##多个AS时得到最保守的AS
# A2 vs At
cat 222|awk '$1~/^evm/&&$2~/^Gh/{print $0}$1~/^Gh/&&$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_At_conserve_end
#D5 vs Dt
cat 222|awk '$1~/^Gor/&&$2~/^Gh/{print $0}$1~/^Gh/&&$2~/^Gor/{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 >D5_Dt_conserve_end
# A2 vs D5
cat 222|awk '$1~/^evm/&&$2~/^G/{print $0}$1~/^Go/&&$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_D5_conserve_end
#At vs Dt
cat 222|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 >At_Dt_conserve_end同源基因间AS event数目统计
统计剪切事件在同源基因间发生的次数
#A2 vs At
python ~/scripte/Alternative/module/homologASeventCount.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt -AS1 ~/work/Alternative/result/homologo/A2/end_third -AS2 ~/work/Alternative/result/homologo/TM-1/end_third -T IntronR -o A2_At_IR
# D5 vs Dt
python ~/scripte/Alternative/module/homologASeventCount.py -homolog ~/work/Alternative/result/homologo/homologGene/D5_vs_Dt_collinearity.txt -AS1 ~/work/Alternative/result/homologo/D5/end_third -AS2 ~/work/Alternative/result/homologo/TM-1/end_third -T IntronR -o D5_Dt_IR
##最后统一去重
tmp=`mktemp`
for i in `ls ./`; do sort ${i}|uniq >${tmp} && cp ${tmp} ${i}; done只有同源基因都存在剪切事件，才会考虑保守性；如果仅仅只有一个同源基因存在，到时候可以看看对应的亚基因组有没有保守的事件
consevePairegenes目录
for i in `ls ./homologeventStatic`; do awk '$2>0&&$4>0{print $0}' ./homologeventStatic/${i} >./consevePairegenes/${i}; done
##对应的事件数目
for i in `ls ./`
do
printf $i"\t"
awk '{a+=$2+$4}END{print a}' ${i}
done
##保守的基因对和事件数目相同亚基因组比较
比较
基因对数
事件数目
保守基因对
保守事件
IR A2 vs At
4508
31401
3075
11618
D5 vs Dt
5574
37375
3852
13812
ES A2 vs At
388
1285
256
542
724
2428
547
1192
AltA A2 vs At
1201
4239
610
1372
1661
5591
1031
2288
AltD A2 vs At
868
3041
430
948
1192
4000
713
1564
> >
在IR事件上，A基因组中共有4508对同源基因存在31401个IR事件，其中3075对基因表现出保守的IR，保守的IR事件占总共事件的36.99%；D基因组中共有5574对同源基因存在37375个IR事件，其中3852对占据13812件（36.95%）保守事件
在ES事件上，A基因组中388对同源基因有1285个事件，其中256对同源基因271（21.08%)为保守的；
D基因组中724对同源基因有2428个事件，其中547对同源基因1192（49.09%）为保守的
在AltA中，A基因组1201对同源基因有4239个事件，其中610对同源基因1372（32.37%）为保守的AltA
D基因组中1661对同源基因有5591个事件，其中1031对基因2288（40.9%）为保守的AltA
在AltD中，A基因组868对同源基因3041个事件，其中430对同源基因31.17%为保守的AltD
D基因组中1192对同源基因4000个事件，其中713对同源基因1564（39.1%)为保守的AltD
不同亚基因组之间进行比较看看
At vs Dt与A2 vs D5 在IR保留事件上
比较
基因对数
事件数目
保守基因对数
事件数目
IR A2 vs D5
4849
32875
2908
10110
At vs Dt
4737
31889
3219
11360
ES A2 vs D5
446
1521
233
494
At vs Dt
462
1502
297
638
AltA A2 vs D5
1136
3840
518
1122
At vs Dt
1527
5217
911
2036
AltD A2 vs D5
729
2409
317
682
At vs Dt
1220
4256
722
1632
IR 30.75%~35.62%
ES: 32.48%~42.48%
AltA:29.22%~39.03%
AltD:28.31%~38.35%
比较A2 D5在二倍体时期，保守的IR，与At、Dt四倍体时期保守的IR看看存在多少交集
使用画图脚本conserveAS.R
保守事件与非保守事件、总事件之间长度比较
python ~/scripte/Alternative/module/extractAsEventCoordiate.py
-p ~/work/Alternative/result/Ga_result/CO11_12_result/06_Alignment/all.collapsed.gtf
-r ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gtf
-AS ~/work/Alternative/result/homologo/A2/end_third
-IR IR_length -ES ES_length -AltA AltA_length -AltD AltD_length
## 提取保守的与不保守的事件的长度
for i in ES IR AltA AltD; do
cat ../blastdb/${i}/A2_At_conserve_end ../blastdb/${i}/A2_D5_conserve_end |cut -f1|sort |uniq |awk '{print "-\t-\t-\t"$0}'|cat - ../A2/${i}_length|sort -k4|uniq -f3 -d >${i}/conserve_${i}.txt;
cat ../blastdb/${i}/A2_At_conserve_end ../blastdb/${i}/A2_D5_conserve_end |cut -f1|sort |uniq |awk '{print "-\t-\t-\t"$0}'|cat - ../A2/${i}_length|sort -k4|uniq -f3 -u >${i}/noconserve_${i}.txt;
done
## 绘图数据
for i in ES AltA AltD IR; do for j in conserve noconserve; do awk '{print "'$i'""_""'$j'""\t"$3-$2+1"\t""'$j'"}' ${i}/${j}_${i}.txt >>./plotData ; done; done
可以看到在AltA、AltD、IR中保守的剪切事件相比与不保守的剪切事件更长；这可能是由于在植物中剪切复合体通过识别intron，对前体mRNA进行剪切，因此intron越长不容易被剪切复合体识别，从在mRNA中被保留；但是由于保守的AS事件，可能存在一定的功能，受到选择压的作用，即使长度相比与不保守的AS仍旧不能够被剪切。
但是保守的ES与不保守的ES在长度上没有显著性的差别。
可参考 Genome-wide survey of Alternative Splicing in Sorghum Bicolor
保守事件中碱基数为3的比例
awk '{a=$3-$2+1;if(a%3==0){b+=1};}END{print b/NR}' conserve_IR.txt事件
保守事件
不保守事件
IR
0.347536
0.335992
AltA
0.372501
0.340563
AltD
0.392737
0.334963
ES
0.408131
0.420226
基因保守率计算
（保守的事件数/A基因组事件数+保守事件数/B基因组事件数）/2
## 计算每对基因保守程度
python ~/scripte/Alternative/module/conservePrecentage.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt -isform ../../homologeventStatic/A2_At_IR -AS A2_At_conserve_end -o conserve_precentage.txt
## 转换格式
awk 'NR%2!=0{printf $1"\t"$2"\t"}NR%2==0{print $1"\t"$2}' conserve_precentage.txt
## 平均保守率
cat conserve_precentage.txt|awk 'NR%2!=0{print $0}' |awk '{a+=$2}END{print a/NR}'
## A2_vs_D5 与At_vs_Dt中同源基因在AS保守性上的比较
0.533165 0.574669 也没有多大区别存在保守剪切事件的基因中，ES事件的平均保守率为0.7856，AltA事件的平均保守率为0.7646，AltD事件的平均保守率为0.7636，IR事件的平均保守率为0.5841。
保守率的高低可能和对应事件碱基数是否是3的倍数有关
提取保守基因对对应的AS事件，在碱基为3的倍数的比例
## 转换一个保守率的格式
awk 'NR%2!=0{printf $1"\t"$2"\t"}NR%2==0{print $1"\t"$2}' AltA
## 运行脚本
python ~/scripte/Alternative/module/BaseConserve.py -length1 ../../A2/AltA_length -length2 ../../D5/AltA_length At_length -length4 ../../TM1/AltA_Dt_length -AsCount1 ../../homologeventStatic/A2_At_AltA -AsCount2 ../../homologeventStatic/A2_D5_AltA -AsCount3 ../../homologeventStatic/At_Dt_AltA -AsCount4 ../../homologeventStatic/D5_Dt_AltA -precentage ./111 -o 222
## 批量脚本
for i in AltA AltD ES IR; do awk 'NR%2!=0{printf $1"\t"$2"\t"}NR%2==0{print $1"\t"$2}' ${i} >111; python ~/scripte/Alternative/module/BaseConserve.py -length1 ../../A2/${i}_length -length2 ../../D5/${i}_length -length3 ../../TM1/${i}_At_length -length4 ../../TM1/${i}_Dt_length -AsCount1 ../../homologeventStatic/A2_At_${i} -AsCount2 ../../homologeventStatic/A2_D5_${i} -AsCount3 ../../homologeventStatic/At_Dt_${i} -AsCount4 ../../homologeventStatic/D5_Dt_${i} -precentage ./111 -o 222; paste ./111 ./222 >./Baseconserve/${i}; done
## 事件碱基数为3的比率和保守率
事件类型 事件保守率 碱基数为3的比率
AltA 0.764633 0.370251
AltD 0.763523 0.376887
ES 0.785635 0.376639
IR 0.584109 0.341613
可能是由于IR事件比较多的原因，所以只用那些一个基因里所有事件都保守的去看保守事件在基因组的位置，所有的剪切事件在基因组的位置
提取gene坐标，制作bin,固定100个窗口
## 制作bin
~/software/bedtools2-2.29.0/bin/windowMaker -i winnum -n 100 -b A2_gene.bed >A2_window.bed
## 提取保守的剪切事件bed文件
cat ../../blastdb/IR/A2_At_conserve_end ../../blastdb/IR/A2_D5_conserve_end|cut -f1|sort |uniq |awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"
## 取交集，统计每个bin所在总的百分比
awk '{if($5=="."){a[$4]+=0}else{a[$4]+=1}}END{for(i in a){print i"\t"a[i]}}' 1 >2看文章好像gene是按照intron来分bin的
做了一下，确实是这样的。
## 提取转录本1的量化坐标
## 取交集
~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ../genomeAnnotion/A2.bed -b ./A2_conserve.bed >1
awk '{if($6=="."){a[$4]+=0}else{a[$4]+=1}}END{for(i in a){print i"\t"a[i]}}' 1 >2
total=`wc -l A2_conserve.bed|awk '{print $1}'`
## 做出50个bin，并且计算占的比例
awk '{for(i=0;i<1;i+=0.02){if($1>=i&&$1<i+0.02){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total'"}}' 2 >3
## 20个bin
awk '{for(i=0;i<1;i+=0.05){if($1>=i&&$1<i+0.05){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total'"}}' 2 >3
## 不仅仅看保守的事件，还有整体趋势
cat ~/work/Alternative/result/homologo/A2/end_third |awk '$3~/IntronR/{print $1,$4,$5}' OFS="\t" |sort -k1,1 -k2,3n|uniq >A2_all.bed
## for 循环
for i in A2 D5 TM1; do ~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ../genomeAnnotion/${i}_isform1.bed -b ./${i}_conserve.bed >1; ~/software/bedtools2-2.29.0/bin/intersectBed -loj -a ../genomeAnnotion/${i}_isform1.bed -b ./${i}_all.bed >2; awk '{if($6=="."){a[$4]+=0}else{a[$4]+=1}}END{for(i in a){print i"\t"a[i]}}' 1 >3; awk '{if($6=="."){a[$4]+=0}else{a[$4]+=1}}END{for(i in a){print i"\t"a[i]}}' 2 >4; total1=`wc -l ${i}_conserve.bed|awk '{print $1}'`; total2=`wc -l ${i}_all.bed|awk '{print $1}'`; awk '{for(i=0;i<1;i+=0.02){if($1>=i&&$1<i+0.02){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total1'"}}' 3 >${i}_conservePosition50.txt; awk '{for(i=0;i<1;i+=0.05){if($1>=i&&$1<i+0.05){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total1'"}}' 3 >${i}_conservePosition20.txt; awk '{for(i=0;i<1;i+=0.02){if($1>=i&&$1<i+0.02){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total2'"}}' 4 >${i}_allPosition50.txt; awk '{for(i=0;i<1;i+=0.05){if($1>=i&&$1<i+0.05){a[i]+=$2}}}END{for(i in a){print i"\t"a[i]"\t"a[i]/"'$total2'"}}' 4 >${i}_allPosition20.txt; done
## 合并画图数据
for i in A2 D5 TM1; do for j in AltA AltD ES IR; do awk '{print $0"\tconserve\t""'${j}'"}' ../${j}/${i}_conservePosition20.txt >>${i}; awk '{print $0"\tall\t""'${j}'"}' ../${j}/${i}_allPosition20.txt >>${i}; done; doneLast updated
Was this helpful?