保守的AS
Copy 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 这个问题之后再看看
Copy ###这种情况得调整一下，两个IR事件，坐标实际上只相差了一个碱基；感觉可以合并后整成一个
- - Gorai.001G041400-Chr01-3801948-3802062 Ghir_D07G004050-Ghir_D07-4143773-4143886
- - Gorai.001G041400-Chr01-3801948-3802062 Ghir_D07G004050-Ghir_D07-4143773-4143887 保守的IR
保守ES
保守AltA
保守AltD
计算IR和ES事件的PIR值
保守事件对应的长度
TE insertion in retained introns may be a significant phenomenon in higher plants.
Copy ## 提取保守的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]的，这类可能是基因没有表达导致的
Copy ##没有表达的基因，从文件内过滤掉
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值的关系
模拟临界PIR值
从所有的intron中抽取1000个intron，用PIR阈值进行判断，得到200个为IR，而其中200个判断的IR有多少比例的是真正的IR，每个PIR处进行1000次迭代
随着PIR阈值的不断增加，真实的IR比例占潜在的比例反而越来越小；
PRI值增加，表明这个conserve AS多倍化后存在功能,conserve_isform_conserve_IR
Copy 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
找到所有同源基因在AS处对应的坐标
Copy ## 这种结果需要过滤一下，理论上是保守的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
Copy ## 在鉴定的保守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 GO功能富集分析
将这些存在保守IR事件的四倍体基因提出来，进行GO富集
Copy ## At和Dt
## A2和At
## D5和Dt 特异性的AS在另外一个基因组的坐标
Copy ##提取一个基因组的所有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
Copy 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的原因
表观上的差异
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；这有可能是染色体交叉互换的结果
