重新开始鉴定AS.md
之前是分别使用对应的三代测序数据比对到各自的参考基因组，现在考虑用三代的测序数据分别比对对各自的亚基因组上。
## 将四倍体基因组分层At和DT进行比对，由于Scaffold没法区分A、D就不要了
grep Ghir_A -A1 Ghirsutum_genome_HAU_v1.0_another.fasta > At.fasta
grep Ghir_D -A1 Ghirsutum_genome_HAU_v1.0_another.fasta > Dt.fasta
## gmap 构建索引
/public/home/zpliu/software/gmap-2019-09-12/bin/gmap_build -D /public/home/zpliu/work/Alternative/data/gmap_build -d Ghirsutum_Dt /public/home/zpliu/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Dt.fasta &
## gmap进行比对
cat ./all.collapsed.rep.fa|~/software/gmap-2019-09-12/bin/gmap -D ~/work/Alternative/data/gmap_build/Ghirsutum_At/ -d Ghirsutum_At -f 3 -t 10 -n 1 --min-trimmed-coverage=0.85 --min-identity=0.85 -p 1 -e >A2_At.gtf
## alternative_splice.py进行鉴定
/usr/bin/python ~/software/pipeline-for-isoseq/other/alternative_splice.py -i ./D5_Dt.gff -g ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Ghirsutum_gene_model_motify_Dt.gtf -f ~/work/Alternative/data/Ghirsutum_genome_HAU_v1.0/Dt.fasta -as -ats T -os -t exon -o ./alter_ten -op通过比较发现，如果真的是很保守的事件的话，他们的长度也是很保守的。对于那些长度摆动的先不管，用gmap来比对然后确定在对于同源基因上的位置。
## 将保守的IR事件，提取在对于基因组上的坐标
##不保守的事件在其他基因组的坐标使用gmap来获取对于AltA、AltD进行同样的操作
## 从剪切文件中获取对应的AltA、AltD事件坐标和长度
awk '$3=="AltA"{print $0}' end_third|paste AltA_stand -|awk '{print $2"\t"$5"\t"$6"\t"$1"\t"$3"\t-\t-\t-\t"$6-$5+1}' >AltA_static.txt
awk '$3=="AltD"{print $0}' end_third|paste AltD_stand -|awk '{print $2"\t"$5"\t"$6"\t"$1"\t"$3"\t-\t-\t-\t"$6-$5+1}' >AltD_static.txt对各种剪切事件在进化上进行分类
## 全都保守的
python ../../script_category/AS_all_converse.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold all_converse_AltA.txt
## 祖先和后代保守
python ../../script_category/As_At_A2_converse.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold At_A2_converse_AltA.txt
python ../../script_category/AS_Dt_D5_converse.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold Dt_D5_converse_AltA.txt
## 祖先与祖先保守
python ../../script_category/AS_A2_D5_converse.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold A2_D5_converse_AltA.txt
python ../../script_category/AS_At_Dt_converse.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold At_Dt_converse_AltA.txt
## 混合的
python ../../script_category/As_mix.py ../homoloe_AltA.txt ../../../GhDt_Gr_GhAt_Ga_end_noScaffold mix_AltA.txt比较保守事件与非保守事件的长度差异
## 提取位置坐标
cat *converse*|awk -F "\t" '$2!="-"{print $2,$3,$4"\n"$9,$10,$11"\n"$16,$17,$18"\n"$23,$24,$25}' OFS="\t"|awk '$1~/[A-Za-z]/{print $0"\t"$3-$2+1}'|sort |uniq >converse.bed
## 提取不保守的AS长度
cut -f1,2,3 homoloe_AltD.txt |sort |uniq|awk '{print $0"\t"$3-$2+1}'|cat - evolution/converse_AltD.bed |sort |uniq -u >noConverse_AltD.bed
## 计算平均值
awk '{a+=$4}END{print a/事件数目}' evolution/converse_AltD.bed事件类型 | 保守类型 | 非保守类型 |
AltD | 365 | 536 |
AltA | 343 | 517 |
ExonS | 116 | 1847 |
IR | 264 | 341 |
比较序列是否存在差异
之前基因文件的版本用错了，然后提出来的序列有些问题，然后今天用新的版本的跑了一下才是一样的；还好我之前跑甲基化数据的时候用对了基因组。
对可变剪切进行分类还是使用剪切位点上下游各300bp组成的序列进行blast
提取剪切位点上下游各300bp的序列，叫做splice anchor sequence tags
## 提取上游300bp序列
awk '{print $1,$2-301,$2-1,$5"_"$1"_"$2"_"$3"_"$9"_"$4}' OFS="\t" ../A2_intronR.txt |sort -k2,3 -n |uniq >left.bed
## 提取下游300bp序列
awk '{print $1,$3,$3+300,$5"_"$1"_"$2"_"$3"_"$9"_"$4}' OFS="\t" ../A2_intronR.txt |sort -k2,3 -n |uniq >right.bed
## 使用bedtools提取序列
~/software/bedtools2-2.29.0/bin/fastaFromBed -fi ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta -fo right.fasta -name -bed right.bed
##合并两端序列
python merge.py left.fasta right.fasta D5_intron_junction.fasta进行blast比对
## 将四个同源基因序列构成一个库,id太长了把 -parse_seqids 参数去掉
makeblastdb -in intron_junction.fasta -dbtype "nucl" -parse_seqids -out intron_junction
## 分别对各个基因组进行比对
blastn -query Dt_intron_junction.fasta -db blastDB/intron_junction -evalue 1e-5 -outfmt 6 -out 111
## 去提取在各个基因组都比较保守的intron事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]>=1&&a[i][2]>=1&&a[i][1]>=1)print i}}'
## 查看某一个事件是不是对的
grep Ghir_D10G003010_2517095_2517293_199_+ Dt_query_blast.txt |awk '$11==0&&$2!~/Ghir_D.*/{print $0}'提取只在某两个基因组上保守的剪切事件
## 提取只在Dt与D5中保守的事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][2]==0&&a[i][1]>=1)print i}}' >Dt_D5_converse_intronR
## 只在Dt与At中保守的剪切事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][2]==0&&a[i][0]>=1)print i}}'
## Dt与A2中保守的剪切事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt|awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][2]>=1&&a[i][0]==0)print i}}'
## A2与At中保守的剪切事件
awk '$11==0&&$2!~/evm.*/{print $0}' A2_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][0]>=1&&a[i][2]==0)print i}}'
## At与D5中的保守剪切事件
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]==0&&a[i][0]==0&&a[i][1]>=1)print i}}'
## A2与D5中的保守剪切事件
awk '$11==0&&$2!~/evm.*/{print $0}' A2_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][1]>=1&&a[i][2]==0)print i}}'在三个基因组上保守的剪切事件
## Dt、D5、A2中都出现
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][2]>=1&&a[i][1]>=1)print i}}'
## Dt、D5、At中都出现
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]>=1&&a[i][2]==0&&a[i][0]>=1)print i}}'
## At、D5、A2中都出现
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]==0&&a[i][0]>=1&&a[i][1]>=1)print i}}'
## At、Dt、A2中都出现
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]>=1&&a[i][0]>=1&&a[i][1]==0)print i}}'只在某一个基因组中出现的剪切事件
## 只在Dt中出现的剪切事件
cat IR_bed/*Dt*|awk -F "\t" '$5~/Ghir_D/{print $5}' |sort |uniq |awk '{print ">"$0}'|cat - Dt_intron_junction.fasta |grep ">"|sort|uniq -u|sed 's/>//g' |awk -F "_" '{print $3"_"$4,$5,$6,$8,$0}' OFS="\t" >noconverse/Dt_IR.bed
## 只在D5中出现的
cat IR_bed/*D5* |awk -F "\t" '$5~/Gor/{print $5}' |sort |uniq |awk '{print ">"$0}'|cat - D5_intron_junction.fasta |grep ">"|sort|uniq -u|sed 's/>//g' |awk -F "_" '{print $2,$3,$4,$6,$0}' OFS="\t" >noconverse/D5_IR.bed
## 只在A2中出现,非保守的剪切事件
cat IR_bed/*A2*|awk -F "\t" '$5~/evm/{print $5}' |sort |uniq |awk '{print ">"$0}'|cat - A2_intron_junction.fasta |grep ">"|sort|uniq -u|sed 's/>//g' |awk -F "_" '{print $2,$3,$4,$6,$0}' OFS="\t" >noconverse/A2_IR.bed
## 只在At中出现
cat IR_bed/*At*|awk -F "\t" '$5~/Ghir_A/{print $5}' |sort |uniq |awk '{print ">"$0}'|cat - At_intron_junction.fasta |grep ">"|sort|uniq -u|sed 's/>//g' |awk -F "_" '{print $3"_"$4,$5,$6,$8,$0}' OFS="\t" >noconverse/At_IR.bed将对于的剪切事件做成bed文件
## 在所有基因组中都保守的IR对应的bed
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]>=1&&a[i][2]>=1&&a[i][1]>=1)print i}}'|xargs -I {} grep {} Dt_query_blast.txt| awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/all_converse.bed
## 提取只在Dt与D5中保守的事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][2]==0&&a[i][1]>=1)print i}}'|xargs -I {} grep {} Dt_query_blast.txt|awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/Dt_D5_converse.bed
## 只在Dt与At中保守的剪切事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][2]==0&&a[i][0]>=1)print i}}' |xargs -I {} grep {} Dt_query_blast.txt|awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/Dt_At_converse.bed
## Dt与A2中保守的剪切事件
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt|awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][2]>=1&&a[i][0]==0)print i}}' |xargs -I {} grep {} Dt_query_blast.txt|awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/Dt_A2_converse.bed
## A2与At中保守的剪切事件
awk '$11==0&&$2!~/evm.*/{print $0}' A2_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]==0&&a[i][0]>=1&&a[i][2]==0)print i}}'|xargs -I {} grep {} A2_query_blast.txt |awk '$11==0&&$2!~/evm.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/A2_At_converse.bed
## At与D5中的保守剪切事件
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]==0&&a[i][0]==0&&a[i][1]>=1)print i}}'|xargs -I {} grep {} At_query_blast.txt|awk '$11==0&&$2!~/Ghir_A.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/At_D5_converse.bed
## A2与D5中的保守剪切事件
awk '$11==0&&$2!~/evm.*/{print $0}' A2_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][1]>=1&&a[i][2]==0)print i}}'|xargs -I {} grep {} A2_query_blast.txt |awk '$11==0&&$2!~/evm.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/A2_D5_converse.bed
## Dt、D5、A2中都出现
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]==0&&a[i][2]>=1&&a[i][1]>=1)print i}}' |xargs -I {} grep {} Dt_query_blast.txt|awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/Dt_D5_A2_converse.bed
## Dt、D5、At中都出现
awk '$11==0&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][1]>=1&&a[i][2]==0&&a[i][0]>=1)print i}}' |xargs -I {} grep {} Dt_query_blast.txt|awk '$11==0&&$2!~/Ghir_D.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/Dt_At_D5_converse.bed
## At、D5、A2中都出现
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]==0&&a[i][0]>=1&&a[i][1]>=1)print i}}'|xargs -I {} grep {} At_query_blast.txt|awk '$11==0&&$2!~/Ghir_A.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/At_D5_A2_converse.bed
## At、Dt、A2中都出现
awk '$11==0&&$2!~/Ghir_A.*/{print $0}' At_query_blast.txt |awk '$2~/evm/{a[$1][0]+=1}$2~/Ghir_D/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][2]>=1&&a[i][0]>=1&&a[i][1]==0)print i}}'|xargs -I {} grep {} At_query_blast.txt|awk '$11==0&&$2!~/Ghir_A.*/{print $1"\n"$2}'|awk -F "_" '$1~/Ghir/{print $3"_"$4,$5,$6,$8,$0}$1!~/Ghir/{print $2,$3,$4,$6,$0}' OFS="\t" |sort |uniq >IR_bed/At_Dt_A2_converse.bed查看基因对的数目
cut -f5 all_converse.bed |awk -F "_" '$1~/Ghir/{print $1"_"$2}$1!~/Ghir/{print $1}' |sort|uniq |grep "Gor" |wc -l事件类型 | 基因对的数目 | 事件数目 |
四个都保守的 | 843 | 1390 |
At = A2 | 1515 | 2629 |
Dt =D5 | 1459 | 2189 |
At=Dt | 1044 | 1528 |
A2=D5 | 1060 | 1511 |
只在At不存在 | 608 | 845 |
只在Dt不存在 | 653 | 945 |
只在A2不存在 | 837 | 1199 |
只在D5不存在 | 574 | 830 |
只有At有 | 4346 | 9387 |
只有Dt有 | 4429 | 9780 |
只有A2有 | 4541 | 10234 |
只有D5有 | 55007 | 12117 |
统计保守的IR长度信息
在统计长度的时候发现 Gorai.006G213800这个基因内含子长度都6000多了，对应的At与A2的内含子长度只有588多但是仍旧都可以发生IR
cat converse/* |sort|uniq |awk -F "_" '{print $(NF-1)}'|less
cat noconverse/* |sort|uniq |awk -F "_" '{print $(NF-1)}'|less219 12 30
师兄说不能只根据 evalue 来筛选保守的片段，得看相似度和覆盖度
## 将相似度与覆盖度的指标写进结果内
blastn -query D5_intron_junction.fasta -db blastDB/intron_junction -outfmt "6 qseqid sseqid qstart qend sstart send nident pident qcovs evalue bitscore" -evalue 1e-5 -out D5_query_blast.txt
## 这个筛选的指标感觉还可以
awk '$7>200&&$8>90&&$9>85{print $0}' At_query.blast.txt把筛选的命令重新写一遍
## 都保守的事件，根据相似度和覆盖度筛选
awk '$7>200&&$8>90&&$9>85&&$2!~/Ghir_D.*/{print $0}' Dt_query_blast.txt |awk '$2~/Ghir_A/{a[$1][0]+=1}$2~/evm/{a[$1][2]+=1}$2~/Gor/{a[$1][1]+=1}END{for(i in a){if(a[i][0]>=1&&a[i][2]>=1&&a[i][1]>=1)print i}}'Last updated