保守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

AS保守与非保守

可以看到在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; done