分析同源基因间可变剪切的差异

筛选各个棉种中发生可变剪切的同源基因

比较剪切位点处甲基化水平

使用单个碱基甲基化程度百分比的程度来衡量甲基化水平

# 使用之前的脚本，把不是0的之间换成1就可以了
对标准化的文件做处理    exon_5scale_single_base.txt
awk '$2==0{print $0}$2!=0{print $1"\t"1}' exon_5scale_single_base.txt
# 再计算每个位点的甲基化概率
for i in 1
do
awk '$2==0{print $0}$2!=0{print $1"\t"1}' exon_5scale_single_base.txt|sed 's/_/\t/g'|awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "ConExon\t""5_exon\t"i"\t"a[i]/87456}' >exon_5_end
awk '$2==0{print $0}$2!=0{print $1"\t"1}' exon_3scale_single_base.txt |sed 's/_/\t/g'|awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "ConExon\t""3_exon\t"i"\t"a[i]/87456}' >exon_3_end
awk '$2==0{print $0}$2!=0{print $1"\t"1}' intron_3scale_single_base.txt |sed 's/_/\t/g'|awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "ConExon\t""3_intron\t"i"\t"a[i]/87456}' >intron_3_end
awk '$2==0{print $0}$2!=0{print $1"\t"1}' intron_5scale_single_base.txt |sed 's/_/\t/g'|awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "ConExon\t""5_intron\t"i"\t"a[i]/87456}' >intron_5_end
done
## 更新坐标轴
awk '{print $1"\t"$2"\t"$3+200"\t"$4}' exon_5_end |sort -k1 -n|cat - intron_5_end |sort -k3,3 -n  >intron_exon_5.txt 
## 将3'内含子延长75 然后合并
awk '{print $1"\t"$2"\t"$3+75"\t"$4}' intron_3_end |cat - exon_3_end |sort -k3 -n >intron_exon_3.txt

### 3'端数据
ggplot(data=data1,aes(x=data1$V3,y=data1$V4,fill=data1$V1,color=data1$V1))+geom_line(size=1)+
    ylim(c(0,0.008))+
    scale_x_continuous(breaks = c(0,75,274),labels = c('-75','','+200'))+
    theme(panel.grid=element_blank(),panel.border = element_blank(),panel.background = element_rect(fill = "white"))+
    geom_vline(xintercept = 75,linetype="dashed",color="blue")+
    theme(axis.line=element_line(size=0.8))+xlab("")+ylab("average mCpG/CpG")+scale_color_manual(values = c("#1e90ff","#ff7675"))+guides(color=guide_legend(title=NULL))
### 5’端数据
ggplot(data=data2,aes(x=data2$V3,y=data2$V4,color=data2$V1))+geom_line(size=1)+ylim(c(0,0.008))+
    scale_x_continuous(breaks = c(0,200,274),labels = c(-200,'','+75'))+
    theme(panel.grid=element_blank(),panel.border = element_blank(),panel.background = element_rect(fill = "white"))+
    geom_vline(xintercept = 200,linetype="dashed",color="blue")+
    theme(axis.line=element_line(size=0.8))+xlab("")+ylab("average mCpG/CpG")+scale_color_manual(values = c("#1e90ff","#ff7675"))+guides(color=FALSE)

A2_8 vs D5_6

棉种	gene 数目	事件数目
TM1	21405	87934
A2	10861	42044
D5	12684	51616

棉种	ES(ExonS)	IR(IntronR)	AltD	AltA	AltP	Other
TM1	4139/3440	53689/18239	10223/7111	11951/8180	5088/3441	2823/2046
A2	2276/1560	28078/9203	4135/2810	5199/3376	1823/1296	1435/949
D5	4280/2958	31756/10581	4900/3427	7055/4731	1819/1303	1805/1209

事件还得重新统计一次

# 获取A2发生了剪切事件的同源基因编号
cut -f2  ~/work/Alternative/result/Ga_result/CO11_12_result/06_Alignment/alter_eight/Altenative_category/end_third|sed /^g/d|sort |uniq >spliceGeneId.txt
cat spliceGeneId.txt homologTD.txt |sort|uniq -d >homologSpliceGeneId.txt
# 根据发生可变剪切的同源基因，找对应的其他基因组的同源基因
cat homologSpliceGeneId.txt |xargs -I {} grep {} ../GhDt_Gr_GhAt_Ga_end_noScaffold >111
mv 111 homologSpliceGeneId.txt

统计有多少同源基因同时发生了可变剪切

# A2中发生可变剪切，在D5中是否也同样发生可变剪切
cat homologSpliceGeneId.txt|xargs -I {} grep {} GhDt_Gr_GhAt_Ga_end_noScaffold >A2_D5_conSplice.txt

分析基因区域甲基化的程度与发生可变剪切的位置之间的关系

## 对每个基因的区域提取200bp的窗口，作标准
cut -f4 A2_D5_At_Dt_conSplice.txt |xargs -I {} grep  {} ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gff |awk '$3~/gene/{print $0}' |awk '{print $1,$4,$5,$7,$9}' OFS="\t"|sed 's/\;.*//g'|awk '{print $1,$2,$3,$4,substr($5,4)}' OFS="\t"  |awk '{for(i=$2;i<=$3;i+=200){print $1"\t"i"\t"i+200"\t"$4"\t"$5"_"(i-$2)/200}}'>A2_gene_location

### 统计甲基化频率和可变剪切频率
~/software/bedtools2-2.29.0/bin/intersectBed -a A2_gene_location -b  ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -loj >A2_gene_Meth.txt
### 提取剪切bed文件
awk '{print $1,$4,$5,$2}' OFS="\t"   ~/work/Alternative/result/Ga_result/CO11_12_result/06_Alignment/alter_eight/Altenative_category/end_third |sed /^c/d >splice.bed 
### 因为基因长度不一样导致窗口数目是不一样的，在求均值的时候不能直接除以基因数目
cut -f5 A2_gene_location |cut -f2 -d _|sort -n |uniq -c >1 # 获取每个窗口的基因数
awk '$1==$6{print $5"\t1"}$1!=$6{print $5"\t0"}' A2_gene_Meth.txt |awk -F "_" '{print $2}' |awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}'|sort -k1,1 -n >2 #获取每个窗口甲基化的数目包括0个甲基化的位置
paste 2 1|awk '{print "MethlationRate\t"$1"\t"$2/$3}' >Methlation.txt
awk '$1==$6{print $5"\t1"}$1!=$6{print $5"\t0"}'  A2_gene_Splice.txt  |awk -F "_" '{print $2}' |awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}'|sort -k1,1 -n|paste - 1|awk '{print "AlterSpliceRate\t"$1"\t"$2/$3}' >AlterSplice.txt

​

将每个基因的长度标准化成100段或者200段，再比较每个区域甲基化和剪切程度

# 将每个基因区域平均分为100份,一个有小数点部分，一个没有小数点还得去掉，不然Bedtools不会处理
cut -f4 A2_D5_At_Dt_conSplice.txt |xargs -I {} grep  {} ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gff |awk '$3~/gene/{print $0}' |awk '{print $1,$4,$5,$7,$9}' OFS="\t"|sed 's/\;.*//g'|awk '{print $1,$2,$3,$4,substr($5,4)}' OFS="\t" >2
#######
awk '{for(i=$2;i<=$3;i+=($3-$2)/100){print $1"\t"i"\t"i+($3-$2)/100"\t"$4"\t"$5"_"(i-$2)/(($3-$2)/100)}}' 2  >gene_100.bed
### 修改版，防止小数和科学计数的存在
awk '{for(i=$2;i<=$3;i+=($3-$2)/100){printf("%s\t%d\t%d\t%s\t%s_%s\n",$1,i,i+($3-$2)/100,$4,$5,(i-$2)/(($3-$2)/100)) }}' 2  >gene_100.bed

分析同源基因之间isform数目的关系

获取每个基因的isform数目

## 统计每个基因在叶片组织中表达的isform数目的分布
ed '1d' end_third |awk '{print $2"\t"$6;print $2"\t"$7}'|sort|uniq |awk '$1~/^Ghir_D/{a[$1]+=1}END{for(i in a)print i"\t"a[i]}'|cut -f2|awk '{a[$1]+=1}END{for(i in a)print "Dt\t"i"\t"a[i]}'|less

所有isform长度的分布信息,计算的是cDNA的长度，把转录本的长度信息也追加到文件中

## 获得每个isform的编号
sed '1d' end_third |awk '{print $2"\t"$6;print $2"\t"$7}'|sort|uniq |awk -F "|" '{print $1}'|awk '$2~/^P/{print $1"\t\\\""$2"\\\""}$2~/^[^P]/{print $1"\t"$2".exon"}'
## 将Pacbio的注释信息与参考基因组注释信息合并
cat  G.arboreum.Chr.v1.0.gff all.collapsed.gff >isform.gff
## 提取发生AltSplice的isform注释信息
sed '1d' end_third | awk '{print $2"\t"$6;print $2"\t"$7}' | sort | uniq | awk -F "|" '{print $1}' | awk '$2~/^P/{print $1"\t\\\""$2"\\\""}$2~/^[^P]/{print $1"\t"$2".exon"}' | cut -f2 | xargs -I {} grep {} isform.gff > splice_isform.gff &

比较同源基因发生剪切的片段的长度

使用bedtool提取对应位置的基因序列

~/software/bedtools2-2.29.0/bin/fastaFromBed -fi ~/genome_data/Ghirsutum_genome_HAU_v1.1/Ghirsutum_genome.fasta   -fo Dt.fasta  -name+ -bed Dt_intronR.txt