统计PB支持的AS
提取那些有PacBio支持的AS事件，并且去除了Scaffold，提取剪切事件所发生的转录本；
并且将那些剪切位置相差几个bp的AS合并了
### 对于IR事件发生在有重叠的坐标区域时，把对应的区域进行合并
python arrangeAS.py -AS ~/work/Alternative/result/Gr_result/CO41_42_result/11_AS/end_splice.txt -r ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.1.gene.gtf -p ~/work/Alternative/result/Gr_result/CO41_42_result/07_annotation/merge.gtf -o D5_AS.txt
统计AS的种类和数目
事件数目
基因数目
IR和其他AS的关系
出现IR的地方同时也似乎容易出现AltA、AltD等情况，分析先有IR还是先有AltA、AltD；
获取剪切事件所在的基因组坐标而不是内含子坐标
提取既有IR又有AltA、AltD的基因出来；分析AltA与IR之间是否有着某种联系
## 获取剪切事件的基因组坐标
python findASposition.py
##分析IR与其他剪切事件的坐标之间的位置
#111 剪切事件在基因组的位置
#22 AltA与IR间的距离
#33 同时出现IR与AltA、AtlD的基因
python IR_AltA.py -AS 111 -o 22 -sim 33
##统计多少AltA与AltD与IR相距小于10bp
grep "AltA" 22|awk '{print $1"-"$5"-"$6}'|sort |uniq |wc -l
grep "AltD" 22|awk '{print $1"-"$5"-"$6}'|sort |uniq |wc -l
grep "AltA" 22 |awk '{a[$1"-"$5"-"$6]=0;if($7<=2){a[$1"-"$5"-"$6]+=1}}END{for(i in a){if(a[i]>0){print i}}}'|wc -l
grep "AltD" 22 |awk '{a[$1"-"$5"-"$6]=0;if($7<=2){a[$1"-"$5"-"$6]+=1}}END{for(i in a){if(a[i]>0){print i}}}'|wc -l
在存同时存在IR与AltA、AltD事件时；AltA、AltD与IR之间间隔小于2bp碱基
将近有27%左右的AltA、AltD与IR的事件的坐标是相差一个碱基的差异
鉴定保守的AS
## 统计同源基因间同时存在AS基因数和事件数目
for i in ExonS IntronR AltA AltD; do python ~/scripte/Alternative/module/homologASeventCount.py -homolog ~/work/Alternative/result/homologo/homologGene/A2_vs_At_collinearity.txt -AS1 ../../conserveAS/allAS/A2_AS.txt -AS2 ../../conserveAS/allAS/TM1_AS.txt -T ${i} -o ${i}; tmp=`mktemp`done; sort ${i}|uniq >${tmp}; mv ${tmp} ${i}; done
awk '$2!=0&&$4!=0{print $0}' IntronR|awk '{a+=$2;b+=$4}END{print a,b,NR}'
##统计保守的基因数
cut -f1 -d"-" IR.txt_end |sort|uniq|wc -l
保守的IR
保守的ES
保守的AltA
保守的AltD
uniq IR在另外一个基因组的坐标
这个uniq的IR在另外一个基因存在三种情况：
在另外一个基因组变成exon，已经进化成exon了
##得到uniq的AS
awk '$3~/IntronR/{print $2,$1,$4,$5}' OFS="-" ../../allAS/A2_AS.txt >A2_all_IR.txt
##取uniq
cut -f1 ../IR.txt_end |cat - A2_all_IR.txt|sed 's/-[-+]//g'|sort |uniq -u
提取基因的固定长度的k-mer序列，使用bwa进行比对，获得uniq IR在另外一个基因组上的坐标；考虑到发生uniq IR的内含子在两个基因组中的长度可能会不一样，同时结合两端FEST序列进行一次锚定
## 取出基因组中两个注释文件的所有的Intron坐标，去冗余
python mergeRefPacBioIntron.py -p ~/work/Alternative/result/Gh_result/CO31_32_result/07_annotation/merge.gtf -r ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model.gtf -o 11
## 去冗余和去除scaffold
awk '$4~/Gh/&&$1!~/Sca/{print $0}' 11 |sed 's/\..*//g' |sort -k1,1 -k2,3n|uniq >TM1_AllIntron.txt
## 将IR与对应基因组的所有Intron进行比较，进行wu-Blast
使用k-mer方法进行实验
首先提取与uniq IR一样长的k-mer片段，实验wu-Blast进行比对，筛选必读长度和IR一样的，并且得分最高的
bsub -q smp -n 1 -J D5_A2 -e kmer.err -o kmer.out -R span[hosts=1] "python ../../k-merBWA.py -AS D5_uniq_IR.txt -gff ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gff -fa1 ~/work/Alternative/data/Gr_genome/Graimondii_221_v2.0.fa -fa2 ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta -ho ~/work/Alternative/result/homologo/homologGene/A2_vs_D5_collinearity.txt"
先进行的是k-mer方法，之后又使用两端FEST方法；看两种方法的交集怎么样
由于FEST取的是左右300bp，因此假阳性情况会比较多，尽量与k-mer的数据为准
##获取wu-Blast的最优的结果
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 >At_uniq_Dt.txt
##统计wu-Blast方法中得到保守，在k-mer中同样保守的数目
cut -f1 At_uniq_Dt.txt |xargs -I {} grep {} ./tmp1594714640/end |wc -l
## 合并两组数据，以k-mer的数据为准
sed '/^$/d' end |awk '{print "B\t"$1"\t"$2}' >../end
awk '{print "A\t"$1"\t"$3}' Dt_uniq_D5.txt |cat - end |sort -k2,2|awk '{print $1"\t"$3"\t"$2}'|uniq -f2|awk '{print $3"\t"$2}' >Dt_uniq_D5_end
将保守的与uniq的合并，统计IR在两个基因组的情况
##合并两次检测的保守IR
cut -f1,3 IR.txt_end|sed 's/-[+-]//g' >1
cat uniqA2/A2_uniq_At_end >>1
awk '{print $2"\t"$1}' uniqAt/At_uniq_A2_end >>1
###去个重
sort -k2,2 1|uniq -f1|awk '{print $2"\t"$1}'|sort -k2,2 |uniq -f1|awk '{print $2"\t"$1}' >2
mv 2 1
##统计IR是否发生
python judgeIR.py -all all_AS.txt -i A2_At/1 -o A2_At/2
统计这个区域是外显子、或者内含子
把所有的区域计算一下read的覆盖情况
### At区域的read覆盖情况
cut -f2 A2_At/1 |cat - At_Dt/1 |cut -f1|sort |uniq |awk -F "-" '{print $2"\t"$3"\t"$4"\t"$1}' >PIR_AS/At_Position
## 计算覆盖的read数目
将这段区域，与PacBio注释文件、参考基因组注释文件进行比较；判断它是否一直是以外显子的形式、或者一直以内含子的形式进行转录；还有一种复杂的情况是Ohter类
提取参考基因组和PacBio每个转录本的注释信息；将这段区域与转录本的注释信息进行比较；如果与内含子或者exon的交集，达到这个片段的90%以上，就认为这个区域是被注释为内含子或者外显子；如果在一个转录本中被注释为exon，另外一些被注释成intron；则认为发生了IR事件。最后统计在一个基因中所有转录本的情况；如果都满足注释为内含子或者都瞒住注释为外显子，就认为该区域为组成型内含子或组成型外显子
r1 p1对应2号文件的第一列
## 对保守的区域进行注释
python ../annotion_conservePosition.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 -AS ./2 -o 3333
## 提取每个区域的IR score值
python ../extractPIR.py -all ../all_IR_score.txt -i 3333 -o 4444
## 统计IR score值的变化
grep 'IR\s+intron' -E 4444 |awk '{a+=$5;b+=$6}END{print a/NR"\t"b/NR}'
统计比例
画一个堆积的条形图：
for i in 1
do
grep 'IR\s+exon' -E 4444 |wc -l
grep 'IR\s+intron' -E 4444 |wc -l
grep 'IR\s+other' -E 4444 |wc -l
grep 'IR\s+IR' -E 4444 |wc -l
grep 'exon\s+IR' -E 4444 |wc -l
grep 'intron\s+IR' -E 4444 |wc -l
grep 'other\s+IR' -E 4444 |wc -l
grep 'IR\s+IR' -E 4444 |wc -l
done
A2 vs At
IR保守率：
D5 vs Dt
A2 vs D5
At vs Dt
，大部分变成了intron，少部分变成exon；四倍体中的IR大部分是由二倍体中的intron转变而来的。
并且在A2中有 1049/16545 的IR发生了外显子化； D5中774/15269 的IR发生外显子化。
计算类型PIR值的变化
可以看出，发生外显子化的IR，它的PIR值很高；而发生intron化的内含子PIR值很低；
由于不同基因组在比较的时候，总read数目存在一个差异，统一将TM1再除以1.35倍
PIR值：比对到intron的read/(intron+两端read)
外显子化的IR、与内含子化的IR在PIR值上的差异；多倍化后外显子化的IR，PIR值显著性的增加
##获取每种类型的PIR值的变化
grep -E "y\s+n" 3333 |grep "other"|cut -f1|awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/A2_PIR.txt >test/5
##看PIR平均值的变化
awk '$5+$6==0{a+=0}$5+$6!=0{a+=$6/($5+$6)}END{print a/NR}' 5
比较PIR值的差异：
由exon 变成IR状态，PIR值仍旧保持一个较高的值，说明这个位置不容易被剪切掉
原先是exon，变成IR；仍旧是发生IR占据主导地位；当然也有剪切类型占据主导地位的
IR变成intron的过程，PIR的值不减少，反而有所增加；可能是有些IR没有被鉴定到；或者存在AltA、AltD的情况导致的
由于多倍化物种在比对的时候，有问题因此在对A2和At间进行比较的时候，误差感觉比较大
## 多倍化前后PIR的变化情况
grep "Chr" A2toAt/A2_IR_At_exon.txt |awk '$5+$6==0{a+=0}$5+$6!=0{a+=$6/($5+$6)}END{print a/NR}'
## 制作绘图数据
# IR to exon
grep "Chr" A2_IR_At_exon.txt |awk '$5+$6==0{print "0\tChr\tIR2exon"}$5+$6!=0{print $6/($5+$6)"\tChr\tIR2exon"}' >>plotdata
grep "Ghir" A2_IR_At_exon.txt |awk '$5+$6==0{print "0\tGhir\tIR2exon"}$5+$6!=0{print $6/($5+$6)"\tGhir\tIR2exon"}' >>plotdata
#exon to IR
grep "Chr" A2_exon_At_IR.txt |awk '$5+$6==0{print "0\tChr\texon2IR"}$5+$6!=0{print $6/($5+$6)"\tChr\texon2IR"}' >>plotdata
grep "Ghir" A2_exon_At_IR.txt |awk '$5+$6==0{print "0\tGhir\texon2IR"}$5+$6!=0{print $6/($5+$6)"\tGhir\texon2IR"}' >>plotdata
# IR to intron
grep "Ghir" A2_IR_At_intron.txt |awk '$5+$6==0{print "0\tGhir\tIR2intron"}$5+$6!=0{print $6/($5+$6)"\tGhir\tIR2intron"}' >>plotdata
grep "Chr" A2_IR_At_intron.txt |awk '$5+$6==0{print "0\tChr\tIR2intron"}$5+$6!=0{print $6/($5+$6)"\tChr\tIR2intron"}' >>plotdata
# intron to IR
grep "Chr" A2_intron_At_IR.txt |awk '$5+$6==0{print "0\tChr\tintron2IR"}$5+$6!=0{print $6/($5+$6)"\tChr\tintron2IR"}' >>plotdata
grep "Ghir" A2_intron_At_IR.txt |awk '$5+$6==0{print "0\tGhir\tintron2IR"}$5+$6!=0{print $6/($5+$6)"\tGhir\tintron2IR"}' >>plotdata
在同一个基因组内进行比较
比如在A2中都是IR，为啥有的多倍化之后变成了外显子、内含子、IR、Other
在At中都是IR，不同来源的IR，在剪切效率上同样不同
Intron Ration score: 比对到内含子上的read数/内含子长度；再把read总数给标准化
## 在A2中是IR，在At中变成IR、exon、intron、Other
grep "y\s+n" -E 3333 |grep other|cut -f1 |awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/A2_PIR.txt >A2toAt/A2IRAtOther.txt
grep "n\s+y" -E 3333 |grep other|cut -f2|awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/TM1_PIR.txt >A2toAt/A2OtherAtIR.txt
## 由exon变成IR
grep "n\s+y" -E 3333 |grep exon |cut -f2 |awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/TM1_PIR.txt >A2toAt/A2_exon_At_IR.txt
## 有IR变成exon
grep "y\s+n" -E 3333 |grep exon |cut -f1|awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/A2_PIR.txt >A2toAt/A2_IR_At_exon.txt
## 由IR变成IR
## 由IR变成intron
grep -E "y\s+n" 3333 |grep intron|cut -f1 |awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/A2_PIR.txt >A2toAt/A2_IR_At_intron.txt
## 有intron变成IR
grep -E "n\s+y" 3333 |grep intron|cut -f2|awk -F "-" '{print $2,$3,$4,$1}' OFS="\t"|xargs -I {} grep {} ../PIR_AS/TM1_PIR.txt >A2toAt/A2_intron_At_IR.txt
## 二倍体祖先数据
awk '{print $7"\tIR2exon"}' A2_IR_At_exon.txt >>A2_plot
awk '{print $7"\tIR2intron"}' A2_IR_At_intron.txt >>A2_plot
awk '$1~/Chr/{print $7"\tIR2IR"}' A2_IR_At_IR.txt >>A2_plot
awk '{print $7"\tIR2other"}' A2IRAtOther.txt >>A2_plot
## 四二倍体数据
awk '{print $7"\texon2IR"}' A2_exon_At_IR.txt >>At_plot
awk '{print $7"\tintron2IR"}' A2_intron_At_IR.txt >>At_plot
awk '$1~/Ghir/{print $7"\tIR2IR"}' A2_IR_At_IR.txt >>At_plot
通过对A2中的IR分析发现，大部分的IR在多倍化之后变成了内含子；根据分类发现外显子化的内含子有着最高的IR score值
通过对At中的IR分析发现，大部分的IR在二倍体中同样是IR；同样有些在二倍体中是exon；在四倍体中变成IR；并且这些内含子化的exon在mRNA中仍旧有较高比例的保留
看一下外显子化的IR有多少交集，并且进行GO富集分析
这种外显子化、内含子化，可能起着调控基因表达的作用
A2多倍化为At时，有一个IR发生了外显子化；而D5与Dt始终是外显子；表达量增加
A2多倍化At时，有一个外显子发生了内含子化；而D5与Dt始终是外显子；表达量减少
A2多倍化过程中外显子化的IR：801个，其中有389个IR在D5中是外显子
D5多倍化过程中外显子化的IR: 389个，其中有173个IR在A2中是外显子
## 获取A2到At中外显子化的基因,看一下在D5中的状态
grep "y\s+n" -E A2_At/3333 |grep exon |cut -f1|xargs -I {} grep {} ../A2_D5/3333 >test.txt
grep exon test.txt|wc -l
cut -f5 test.txt|sort|uniq -c
## 获取D5到Dt中外显子化IR，看一下在A2中的状态
grep "y\s+n" -E ../D5_Dt/3333 |grep exon |cut -f1|xargs -I {} grep {} ../A2_D5/3333 >test.txt2
grep exon test.txt2 |wc -l
cut -f5 test.txt2|sort|uniq -c
At基因组中由exon转变为IR：552个，其中有233个在Dt中是exon的状态
Dt基因组中由exon转变为IR：492个，其中有308个在At中是exon的状态
## 获取A2到At过程中，exon转变为IR
grep -E "n\s+y" ../A2_At/3333 |grep exon|cut -f2 |xargs -I {} grep {} 3333 >test.txt
cut -f5 test.txt |sort|uniq -c
## 获取D5到Dt过程中，exon转变为IR
grep "n\s+y" -E ../D5_Dt/3333 |grep exon|cut -f2 |xargs -I {} grep {} 3333 >test.txt2
cut -f5 test.txt2 |sort|uniq -c
GO富集分析
A基因组外显子化的基因、从外显子转变成IR的基因；因为即使变成IR之后PIR值也还是很高的
亚基因组特异性的AS
At与Dt进行比较的时候，例如At中存在IR、Dt中不存在IR；并且At中那个IR的地方有很多read覆盖，而Dt中几乎没有read覆盖
提取At中发生IR的IRscore ，Dt中没有发生IR的IR score；存在显著性差异
At中发生IR、Dt中没有发生IR；并且发生IR的IR score值大于0.1，没有发生IR的值小于0.1
排除那些由于PacBio没检测到的IR
## 筛选亚基因组中特异的IR
grep "intron\s+IR" -E test.txt|awk '$6>$5&&$5<=0.1{print $0}'|wc -l
## 提取所有事件的分类，和对应的IRS值
python extractPIR.py -all all_IRScore.txt -i At_Dt/3333 -o At_Dt/test.txt
