参考： 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范围
Copy ##内部经过筛选除了，选择两端长度大于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
同一种类型的剪切、不同亚基因组之间建库
Copy ## 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数目统计
统计剪切事件在同源基因间发生的次数
Copy #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目录
Copy 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事件上，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 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
保守事件与非保守事件、总事件之间长度比较
Copy 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的比例
Copy awk '{a=$3-$2+1;if(a%3==0){b+=1};}END{print b/NR}' conserve_IR.txt 基因保守率计算
（保守的事件数/A基因组事件数+保守事件数/B基因组事件数）/2
Copy ## 计算每对基因保守程度
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的倍数的比例
Copy ## 转换一个保守率的格式
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个窗口
Copy ## 制作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的
做了一下，确实是这样的。
Copy ## 提取转录本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 
