01全长转录组数据处理
测序数据如下图

每个棉种含有两个重复的数据,之所以有3个文件,是因为有一个重复的数据量没有测够;之后进行了补测
2019-10-10
合并两次测序的bam文件
## 分别对两个bam文件进行排序,按照序列编号顺序排序
samtools sort -n -@ 2 ./A2/dT_BC5_subreads.bam -O BAM -o ./A2/dT_BC5_subreads_sorted.bam
## 将两个排序好的文件进行合并
samtools merge -n -@ 2 -O BAM ./A2/A2_R2.subreads.bam ./A2/dT_BC5_subreads_sorted.bam ./A2/m54139_180604_080709.subreads_sorted.bam
后来发现还得使用ccs方法将两个bam文件指定为一个数据集
2019-10-11
因为在进行数据打磨时报错了
Missing .pbi file
dataset create A2_R2.subreadset.xml --type SubreadSet --generateIndices ./../dT_BC5_subreads.bam ./../m54139_180604_080709.subreads.bam
# 之后使用A2_R2.subreadset.xml文件代替bam文件进行ccs即可
2019-10-13
使用iso-seq中dataset程序处理,重复二的数据
dataset create R2/TM1_R2.subreadset.xml --type SubreadSet --generateIndices R1801371_QJ_BC1_subreads.bam m54139_180607_052119.subreads.bam &
dataset create R2/D5_R2.subreadset.xml --type SubreadSet --generateIndices m54136_180730_021327.subreads.bam m54139_180609_044201.subreads.bam &
使用ccs处理
ccs --noPolish --minPasses 1 --minLength 300 --minSnr 4 --maxLength 10000 --maxDropFraction 0.8 --minPredictedAccuracy 0.8 --numThreads 20 --logFile ./TM1_R2_ccs.log --reportFile ./TM1_R2port_ccs.txt ./TM1_R2.subreadset.xml ./TM1_R2_ccsout.bam
之前数据的分析,师姐已经完成了;接下来主要是针对数据进行个性化的分析了;其实我还挺想自己把流程走一遍的,可惜公司不会把分析流程的protocol发给我们
2019-10-14
将之前的结果解压,到当前文件夹
unrar x result.rar
参考我之前写的软件使用说明
tps://zpliu.gitbook.io/booknote/ke-bian-jian-qie/ruan-jian-shi-yong/01-san-dai-ce-xu-isoseq#2-shi-yong-ccs-dui-yuan-shi-shu-ju-jin-hang-guo-lv
统计测序序列情况
1.二倍体亚洲棉
数据来源:第03三步分析流程里文件Classify_stat.xls
CO_31
序列类型
数目
比例
consensus reads质控后经过ccs处理
320874
100%
non-full-length reads
34265
10.68%
full-length reads(包括chimeric reads)
286277
89.22%
filtered short reads长度小于300bp
332
0.1%
chimeric reads
9710
3.03%
CO_32
序列类型
数目
比例
consensus reads质控后经过ccs处理
292447
100%
non-full-length reads
34265
10.68%
full-length reads(包括chimeric reads)
228041
77.98%
filtered short reads长度小于300bp
366
0.13%
chimeric reads
11003
3.77%
四倍体陆地棉
Col31&C0l32
序列类型
数目
比例
consensus reads质控后经过ccs处理
685383
100%
non-full-length reads
103674
15.13%
full-length reads(包括chimeric reads)
581611
84.86%
filtered short reads长度小于300bp
98
0.01%
chimeric reads
18914
2.76%
二倍体雷蒙德氏棉
Col31&C0l32
序列类型
数目
比例
consensus reads质控后经过ccs处理
487063
100%
non-full-length reads
87450
17.95%
full-length reads(包括chimeric reads)
399112
81.94%
filtered short reads长度小于300bp
501
0.1%
chimeric reads
26652
5.47%
比较测序数据与参考基因组数据中isform的长度与对应exon的数目
在测序数据中
06_Alignment/all.collapsed.gff
文件中有每条isform的信息1.统计每条isform的长度信息,只包含外显子的长度
2.统计每条isform的外显子数目
sed 's/gene_id [^t]*transcript_id//g' all.collapsed.gff |awk '$3~/^e/{isformLen[$9]+=$5-$4+1;exonCount[$9]+=1}END{for(i in isformLen){printf i"\t"isformLen[i]"\t"exonCount[i]"\n"}}'
##得到每个isform的长度和外显子数目 "PB.1024.1"; 4680 5 "PB.10242.1"; 3658 7 "PB.1024.2"; 4463 5 "PB.10243.1"; 2680 3
参考基因组gff文件中有每个mRNA的注释信息
多了一个sed的处理主要是isform的第九列多了·一些信息
ID=evm.model.Ga01G0001.exon2
,把exon那部分去除后获得isform编号awk -F ";" '{print $1}' Ghirsutum_gene_model.gff3|sed 's/\.exon.*//g' |awk '$3~/^e/{isformLen[$9]+=$5-$4+1;exonCount[$9]+=1}END{for(i in isformLen){printf i"\t"isformLen[i]"\t"exonCount[i]"\n"}}'
准备绘图的数据
# 将isform名字信息修改成分组信息 sed -i 's/\"PB[^\t]*/TM1_PBisform/g' my_isform_length_exonCount.txt # 修改基因组的isform信息 sed -i 's/\ID=[^\t]*/TM1_genome/g' my_genome_isform_length_exonCount.txt # 将多个棉种的数据合并成一个文件,比较外显子数目上的差异 awk '{exonClass[$1"*"$3]+=1}END{for(i in exonClass){split(i,tmp,"*");printf tmp[1]"\t"tmp[2]"\t"exonClass[i]"\n"}}' all_genome_isform_comprasion.txt >all_genome_isform_exonCount.txt
得到每个基因组中isform的外显子数目分布情况
#基因组 外显子数目 对应的isform数目 TM1_genome 64 3 D5_PBisform 43 7 D5_genome 16 674 TM1_genome 65 1 D5_PBisform 44 13 D5_genome 30 51
2019-10-16
绘制isform长度比较信息,和外显子数目的比较信息
#绘制isform长度信息 ggplot(data = tmp,aes(x=tmp$V2,y=tmp$V3,fill=tmp$V2))+geom_boxplot()+ facet_wrap(~tmp$V1,scales = "free")+
# 绘制isform中外显子的数目信息
ggplot(data = tmp,aes(x=tmp$V3,y=tmp$V4,fill=tmp$V2))+geom_bar(stat = "identity", position = 'dodge')+scale_x_continuous(breaks = seq(1,20,1))+facet_wrap(~tmp$V1)
2019-10-17
第7个文件夹annotation中包含了,与参考基因组注释文件进行合并后的注释信息
利用地6步分析中得到的非冗余的isform序列,与基因组进行比较
棉种
非冗余的isform
比对到注释基因的isform
没有比对到基因区域
对应基因数
对应的转录本数
原参考注释基因数目
TM1
89,411
83,281
6,130
32,003
39,098
70199
A2
7,2393
69,179
3,214
20,907
参考基因组没有预测不同转录本
40960
D5
55,381
53,058
2,323
18,904
24,368
37505
第8个和第9个、10个文件的数据暂时用不到
分析第11个关于可变剪切的文件
分析不同棉种中可变剪切的类型,及发生可变剪切的基因数目
grep ExonS splice.as.xls|wc -l
统计外显子跳跃事件的数目
grep ExonS splice.as.xls|awk '{print $2}' |sort|uniq |wc -l
统计对应的基因数目
剪切事件数目/对应的基因数目
棉种
gene 总数目
事件总数目
TM-1
23137
121567
A2
11733
60921
D5
14483
81392
棉种
ES(ExonS)
IR(IntronR)
AltD
AltA
AltP
Other
TM1
5824/3519
70132/18813
13147/7329
15631/8522
8182/4201
8650/3621
A2
2536/1564
38121/9515
4746/2736
6129/3278
4764/2318
4624/1975
D5
5502/3005
42372/11062
8104/4713
11799/6094
5980/3046
7634/2810
看看比例TM-1
事件类型
事件比例与数目
基因数目与比例
IR
70132~57%
18813~81%
AD
13147~10.8%
7329~31.6%
AA
15631~12.8%
8522~36.8%
ES
5824~4.8%
3519~15.2%
AP
8182~6.7%
4201~18.1%
Other
8650~7.1%
3621~15.6%
看看比例A2
事件类型
事件比例与数目
基因数目与比例
IR
38121~62.5%
9515~81%
AD
4746~7.8%
2736~23.3%
AA
6129~10%
3278~27.9%
ES
2536~4.1%
1564~13.3%
AP
4764~7.8%
2318~19.7%
Other
4624~7.6%
1975~16.8%
看看比例D5
事件类型
事件比例与数目
基因数目与比例
IR
42372~52%
11062~76%
AD
8104~9.9%
4713~32.5%
AA
11799~14.5%
6094~42%
ES
5502~6.7%
3005~20.7%
AP
5980~7.3%
3046~21%
Other
7634~9.4%
2810~19.4%
D5的剪切事件相对与A2来说有点偏高,可能是由于A2的参考基因组本身就没有对多个转录本进行预测;为了验证这个猜想,就只看PBc得到的转录本之间的剪切事件
棉种
ES(ExonS)
IR(IntronR)
AltD
AltA
AltP
Other
TM1
1225/916
18218/5901
2532/1641
3088/1926
2797/1532
2609/1264
A2
1253/840
17865/5347
2201/1444
2945/1803
2358/1312
2571/1068
D5
504/368
9578/3160
1112/770
1438/901
1379/858
1466/735
# 批量统计次数
for i in 1; do
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep ExonS|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep ExonS|awk '{print $2}' |sort|uniq |wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep IntronR|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep IntronR|awk '{print $2}' |sort|uniq |wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltD|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltD|awk '{print $2}' |sort|uniq |wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltA|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltA|awk '{print $2}' |sort|uniq |wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltP|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep AltP|awk '{print $2}' |sort|uniq |wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep Other|wc -l;
awk '$6~/^PB/&&$7~/^PB/{print $0}' splice.as.xls |grep Other|awk '{print $2}' |sort|uniq |wc -l; done
2019-10-18
在第7个文件夹中matchAnnot_result.txt
文件中包含比对到参考基因组上的isform和对应的参考基因编号
分析全长转录组的数据与原始数据相比在isform数目上的差异
2019-11-2
根据参考基因组的注释文件,产生剪切位点的bed文件
假定剪切复合体是以外显子识别的模式进行剪切的,并且不关心外显子在正负链上的关系
#将第9列改造成只剩下基因ID,结合exon的起始位点,给每个exon一个编号,不分析脚手架了
sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}' >tmp.bed
awk '$1~/Chr/{print $1"\t"$2"\t"$3"\t"$4"\t"$5"_"$2}' constitutiveExon.bed >tmp.bed
### 结果如下
Ghir_A01 13133 13156 - Ghir_A01G000010.1_13133
Ghir_A01 13840 14080 - Ghir_A01G000010.1_13840
Ghir_A01 14164 15594 - Ghir_A01G000010.1_14164
生成各种类型的剪切
# 生成5'外显子区域的bed文件,覆盖到每个碱基 awk '{for(i=$2-1;i<=$2+73;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+1}}' tmp.bed >exon_5.bed awk '{for(i=$3-73;i<=$3+1;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3+73}}' tmp.bed >exon_3.bed awk '{for(i=$2-201;i<=$2-2;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+201}}' tmp.bed >intron_5.bed awk '{for(i=$3+2;i<=$3+201;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3-2}}' tmp.bed >intron_3.bed # 用于提取指定区域的fasta序列 awk '{print $1"\t"$2-1"\t"$2+73"\t"$4"\t"$5}' tmp.bed >exon_5_fasta.bed awk '{print $1"\t"$3-73"\t"$3+1"\t"$4"\t"$5}' tmp.bed >exon_3_fasta.bed awk '{print $1"\t"$2-201"\t"$2-2"\t"$4"\t"$5}' tmp.bed >intron_5_fasta.bed awk '{print $1"\t"$3+2"\t"$3+201"\t"$4"\t"$5}' tmp.bed >intron_3_fasta.bed # for 循环执行 for i in 1 do sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{for(i=$2-1;i<=$2+73;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+1}}' >exon_5.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{for(i=$3-73;i<=$3+1;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3+73}}' >exon_3.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{for(i=$2-201;i<=$2-2;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+201}}' >intron_5.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{for(i=$3+2;i<=$3+201;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3-2}}' >intron_3.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{print $1"\t"$2-1"\t"$2+73"\t"$4"\t"$5}' >exon_5_fasta.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{print $1"\t"$3-73"\t"$3+1"\t"$4"\t"$5}' >exon_3_fasta.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{print $1"\t"$2-201"\t"$2-2"\t"$4"\t"$5}' >intron_5_fasta.bed sed 's/\;.*//g' ../Ghirsutum_gene_model.gtf |sed 's/transcript_id //g' |sed 's/\"//g'|awk '$1~/Gh/&&$3~/exon/{print $1"\t"$4"\t"$5"\t"$7"\t"$9"_"$4}'|awk '{print $1"\t"$3+2"\t"$3+201"\t"$4"\t"$5}' >intron_3_fasta.bed done
将甲基化的数据构造成bed文件
awk '{print $1"\t"$2"\t"$2}' ~/WMJ_3D_DNAmethy/Gh_bismark/CpG_context_D1_binom_fdr_result.txt >CpG_context_D1.bed
统计每个区域中含有甲基化的位点的数目
## 如果区域没有对应的甲基化位点则会显示为-1,有则显示对于的甲基化位置 ~/software/bedtools2-2.29.0/bin/intersectBed -b CpG_context_D1.bed -a exon_3.bed -loj >exon_3Meth_count.txt ## 统计每个区域甲基化的次数,后面几列的坐标一样的话就认为是一个甲基化位点 # awk '$2==$7{print $0}' exon_3Meth_count.txt ## 统计每个区域中单个碱基精度的甲基化程度,需要重新制作bed文件
提取对应区域的CG类型的数目,进行甲基化水平的标准化
# 统计每个区域CG含量水平 python extract_CG_count.py sequencFile outputFile # 提取甲基化是单个碱基的精度下的甲基化水平 python scale_Methlation_by_sigleBase.py exon_3_CGcount.txt exon_3Meth_count.txt scale_single_base.txt # 统计外显子的数目 grep ">" exon_3.fasta|wc -l # 均一化所有外显子的单碱基水平 sed 's/_/\t/g' scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]/外显子数目}' >end
画两幅内含子和外显子连接的图,intron5-exon5 vs exon3-intron3
## 将5'外显子坐标延长200,然后合并 awk '{print $1+200"\t"$2}' exon5_end |sort -k1 -n|cat - intron5_end |sort -k1 -n >intron_exon_5.txt ## 将3'内含子延长75 然后合并 awk '{print $1+75"\t"$2}' intron3_end |cat - exon3_end |sort -k1 -n >intron_exon_3.txt
ggplot2绘制图形
## 5'端 ggplot(data=data2,aes(x=data2$V1,y=data2$V2))+ geom_line(color=c("#ff4757"),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 = "#f1f2f6"))+ geom_vline(xintercept = 200,linetype="dashed",color="blue")+ theme(axis.line=element_line(size=0.8))+xlab("")+ylab("average mCpG/CpG") ## 3'端 ggplot(data=data1,aes(x=data1$V1,y=data1$V2))+ geom_line(color=c("#ff4757"),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 = "#f1f2f6"))+ geom_vline(xintercept = 75,linetype="dashed",color="blue")+ theme(axis.line=element_line(size=0.8))+xlab("")+ylab("average mCpG/CpG")
2019-11-08
重新分析可变剪切的数据,使用gmap进行比对,过滤指标覆盖度85%,相似度90%
cat ./all.collapsed.rep.fa|~/software/gmap-2019-09-12/bin/gmap -D ~/work/Alternative/data/gmap_build/Graimondii_221_v2.0 -d Graimondii_221_v2.0 -f samse -t 10 -n 2 --min-trimmed-coverage=0.85 --min-identity=0.9 >test.sam
2019-11-14
考虑每个isform直接的比较
棉种
ES(ExonS)
IR(IntronR)
AltD
AltA
AltP
Other
TM1
17702/3440
254763/18239
37147/7111
48642/8180
13878/3441
8005/2046
A2
5074/1390
94222/9164
12407/2762
17002/3382
11798/2650
4310/1213
D5
18554/2949
148092/10583
19095/3523
31982/4988
5508/1406
5199/1279
只考虑基因坐标
棉种
ES(ExonS)
IR(IntronR)
AltD
AltA
AltP
Other
TM1
4139/3440
53689/18239
10223/7111
11951/8180
5088/3441
2823/2046
A2
1656/1390
27223/9164
3839/2762
5021/3382
4428/2650
1757/1213
D5
3863/2949
31532/10583
4903/3523
7343/4988
1892/1406
1758/1279
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
2019-11-16
比较Alt exon与constituents exon的甲基化水平
提取可变剪切文件中所有isform的exon坐标
先提取PacBio对应的isform编号,以及对应的mRNA对应的编号
awk '$6~/PB/{print $6"\t"$2}' end_third |sort|uniq >PacBio_ID
awk '$7~/PB/{print $7"\t"$2}' end_third |sort|uniq >>PacBio_ID
awk '$6~/evm/{print $6"\t"$2}' end_third |sort|uniq >mRNA_ID
awk '$7~/evm/{print $7"\t"$2}' end_third |sort|uniq >>mRNA_ID
## 还得去一次重
sort PacBio_ID|uniq >PacBio_ID_1
## 使用python脚本将每个isform的外显子信息坐标信息进去取交集,假定在2个isform中出现的exon坐标,认为是constituents exon
python constitutiveExonLocation.py PacBio_ID mRNA_ID ../test.gff ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gff constitutiveExonLocation
# 给每个constitutive exon取个名字
# 最终的constitutive exon结果里可能会有重复的坐标,由于不同转录本之间太相似了,把一些重叠区域尽量最小化
sort -n -k3,3 constitutiveExon.bed |awk '{print $2,$3,$1,$4,$5}' OFS="\t"|uniq -f 1|awk '{print $2,$1,$3,$4,$5}' OFS="\t" |sort -k2,2 -n |uniq -f1 |awk '{print $3,$2,$1,$4,$5}' OFS="\t"|awk '{print $1,$2,$3,$4,$5"_"$2}' OFS="\t" >tmp.bed
# 提取对应的bed坐标文件,这有点问题没有根据链来提,会有误差
for i in 1
do
awk '{for(i=$2-1;i<=$2+73;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+1}}' tmp.bed >exon_5.bed
awk '{for(i=$3-73;i<=$3+1;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3+73}}' tmp.bed >exon_3.bed
awk '{for(i=$2-201;i<=$2-2;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$2+201}}' tmp.bed >intron_5.bed
awk '{for(i=$3+2;i<=$3+201;i++){print $1"\t"i"\t"i"\t"$4"\t"$5"_"i-$3-2}}' tmp.bed >intron_3.bed
awk '{print $1"\t"$2-1"\t"$2+73"\t"$4"\t"$5}' tmp.bed >exon_5_fasta.bed
awk '{print $1"\t"$3-73"\t"$3+1"\t"$4"\t"$5}' tmp.bed >exon_3_fasta.bed
awk '{print $1"\t"$2-201"\t"$2-2"\t"$4"\t"$5}' tmp.bed >intron_5_fasta.bed
awk '{print $1"\t"$3+2"\t"$3+201"\t"$4"\t"$5}' tmp.bed >intron_3_fasta.bed
done
# 提取对应的序列用于计算CG含量
python ~/scripte/according_CDS_location_find_fasta_mesage.py exon_3_fasta.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta exon_3.fasta
python ~/scripte/extract_CG_count.py exon_3.fasta exon_3_CG_count.txt
# 使用intersectBed做交集
~/software/bedtools2-2.29.0/bin/intersectBed -b CpG_context_D1.bed -a exon_3.bed -loj >exon_3Meth_count.txt
# 提取单碱基下的甲基化水平
python scale_Methlation_by_sigleBase.py exon_3_CGcount.txt exon_3Meth_count.txt scale_single_base.txt
# 按exon的数目进行标准化 # 均一化所有外显子的单碱基水平
grep ">" exon_3.fasta|wc -l
sed 's/_/\t/g' scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "AltExon\t"i"\t"a[i]/外显子数目}' >end
### 搞个for循环懒得敲了
for i in 1
do
python ~/scripte/according_CDS_location_find_fasta_mesage.py exon_3_fasta.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta exon_3.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py exon_5_fasta.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta exon_5.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py intron_3_fasta.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta intron_3.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py intron_5_fasta.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta intron_5.fasta
python ~/scripte/extract_CG_count.py exon_3.fasta exon_3_CG_count.txt
python ~/scripte/extract_CG_count.py exon_5.fasta exon_5_CG_count.txt
python ~/scripte/extract_CG_count.py intron_3.fasta intron_3_CG_count.txt
python ~/scripte/extract_CG_count.py intron_5.fasta intron_5_CG_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -a exon_3.bed -loj >exon_3Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -a exon_5.bed -loj >exon_5Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -a intron_5.bed -loj >intron_5Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -a intron_3.bed -loj >intron_3Meth_count.txt
python ~/scripte/scale_Methlation_by_sigleBase.py exon_3_CG_count.txt exon_3Meth_count.txt exon_3scale_single_base.txt
python ~/scripte/scale_Methlation_by_sigleBase.py exon_5_CG_count.txt exon_5Meth_count.txt exon_5scale_single_base.txt
python ~/scripte/scale_Methlation_by_sigleBase.py intron_5_CG_count.txt intron_5Meth_count.txt intron_5scale_single_base.txt
python ~/scripte/scale_Methlation_by_sigleBase.py intron_3_CG_count.txt intron_3Meth_count.txt intron_3scale_single_base.txt
done
### 最终结果
grep ">" intron_5.fasta |wc -l
## exon的类型记得改一下
for exonCount in 外显子的数目
do
sed 's/_/\t/g' exon_3scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "AltExon\t"i"\t"a[i]/"'$exonCount'"}' >exon_3_end
sed 's/_/\t/g' exon_5scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "AltExon\t"i"\t"a[i]/"'$exonCount'"}' >exon_5_end
sed 's/_/\t/g' intron_5scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "AltExon\t"i"\t"a[i]/"'$exonCount'"}' >intron_5_end
sed 's/_/\t/g' intron_3scale_single_base.txt |awk '{print $(NF-1)"\t"$NF}'|awk '{a[$1]+=$2}END{for(i in a)print "AltExon\t"i"\t"a[i]/"'$exonCount'"}' >intron_3_end
done
### 为了画图方便将所有exon做成连续的坐标
## 将5'外显子坐标延长200,然后合并
awk '{print $1"\t"$2+200"\t"$3}' exon_5_end |sort -k2 -n|cat - intron_5_end |sort -k2 -n >intron_exon_5.txt
## 将3'内含子延长75 然后合并
awk '{print $1"\t"$2+75"\t"$3}' intron_3_end |cat - exon_3_end |sort -k2 -n >intron_exon_3.txt
提取Alter exon的坐标
### 没有染色体正负链的信息,在后面提取序列要用,但是提取CG碱基的时候,由于对称性所以提取哪条链都没关系
### 都以正链为标准
awk '$3~/IntronR/||$3~/ExonS/{print $1,$4,$5,"+",$2}' OFS="\t" end_third >../Alter_exon/alter_exon.bed
2019-11-20
之间比较Constitutive Exon和Alternative Exon之间甲基化水平的差异
Alternative Exon需要把对应的链的信息给补充完整
直接根据基因区编号来补充把
## 提取基因正负链信息
cut -f5 alter_exon.bed |xargs -I {} grep {} ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.gff|awk '$3~/gene/{print $7}' |paste alter_exon.bed - >11111111111
## 过滤重复的可变exon
awk '{print $1,$2,$3,$6,$5}' OFS="\t" ../11111111111 |sort -n -k3,3 |awk '{print $2,$3,$1,$4,$5}' OFS="\t"|uniq -f 1|awk '{print $2,$1,$3,$4,$5}' OFS="\t" |sort -k2,2 -n |uniq -f1 |awk '{print $3,$2,$1,$4,$5}' OFS="\t"|awk '{print $1,$2,$3,$4,$5"_"$2}' OFS="\t" >tmp.bed
## 构造每个外显子在5'|3'端内含子和外显子的甲基化数据
for i in 1
do
awk '$4~/+/{print $1"\t"$2"\t"$2+74"\t"$4"\t"$5}$4~/-/{print $1"\t"$3-74"\t"$3"\t"$4"\t"$5}' ../tmp.bed >exon_5.bed
awk '$4~/+/{print $1"\t"$3-74"\t"$3"\t"$4"\t"$5}$4~/-/{print $1"\t"$2"\t"$2+74"\t"$4"\t"$5}' ../tmp.bed >exon_3.bed
awk '$4~/+/{print $1"\t"$3+1"\t"$3+75"\t"$4"\t"$5}$4~/-/{print $1"\t"$2-75"\t"$2-1"\t"$4"\t"$5}' ../tmp.bed >intron_3.bed
awk '$4~/+/{print $1"\t"$2-75"\t"$2-1"\t"$4"\t"$5}$4~/-/{print $1"\t"$3+1"\t"$3+75"\t"$4"\t"$5}' ../tmp.bed >intron_5.bed
done
## 长度都是一定的,提取对应的CG甲基含量
for i in 1
do
python ~/scripte/according_CDS_location_find_fasta_mesage.py exon_3.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta exon_3.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py exon_5.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta exon_5.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py intron_3.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta intron_3.fasta
python ~/scripte/according_CDS_location_find_fasta_mesage.py intron_5.bed ~/work/Alternative/data/Ga_genome/G.arboreum.Chr.v1.0.fasta intron_5.fasta
done
for i in 1
do
python ~/scripte/extract_CG_count.py exon_3.fasta exon_3_CG_count.txt
python ~/scripte/extract_CG_count.py exon_5.fasta exon_5_CG_count.txt
python ~/scripte/extract_CG_count.py intron_3.fasta intron_3_CG_count.txt
python ~/scripte/extract_CG_count.py intron_5.fasta intron_5_CG_count.txt
done
## 统计每个区域甲基化的碱基数
for i in 1
do
~/software/bedtools2-2.29.0/bin/intersectBed -a exon_3.bed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -loj|awk '$6~/\./{print $5"\t"0}$6~/[^\.]/{print $5"\t"1}'|awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}' >exon_3Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -a exon_5.bed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -loj|awk '$6~/\./{print $5"\t"0}$6~/[^\.]/{print $5"\t"1}'|awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}' >exon_5Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -a intron_3.bed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -loj|awk '$6~/\./{print $5"\t"0}$6~/[^\.]/{print $5"\t"1}'|awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}' >intron_3Meth_count.txt
~/software/bedtools2-2.29.0/bin/intersectBed -a intron_5.bed -b ~/work/Alternative/data/Ga_genome/test/CpG_context_D4.bed -loj|awk '$6~/\./{print $5"\t"0}$6~/[^\.]/{print $5"\t"1}'|awk '{a[$1]+=$2}END{for(i in a)print i"\t"a[i]}' >intron_5Meth_count.txt
done
## 构造ggplot绘图的数据
for i in 1
do
awk '{print "ConExon\t""exon_3\t"$2}' exon_3_CG_count.txt >conExon_CGcount_end
awk '{print "ConExon\t""exon_5\t"$2}' exon_5_CG_count.txt >>conExon_CGcount_end
awk '{print "ConExon\t""intron_3\t"$2}' intron_3_CG_count.txt >>conExon_CGcount_end
awk '{print "ConExon\t""intron_5\t"$2}' intron_5_CG_count.txt >>conExon_CGcount_end
awk '{print "ConExon\t""exon_3\t"$2}' exon_3Meth_count.txt >conExon_MetCount_end
awk '{print "ConExon\t""exon_5\t"$2}' exon_5Meth_count.txt >>conExon_MetCount_end
awk '{print "ConExon\t""intron_5\t"$2}' intron_5Meth_count.txt >>conExon_MetCount_end
awk '{print "ConExon\t""intron_3\t"$2}' intron_3Meth_count.txt >>conExon_MetCount_end
done
感觉失败了

把四个坐标中甲基化位点为0的exon去除试试
#剔除在 exon_5|3 和intron3|5都没有甲基化的exon
cat exon_3Meth_count.txt exon_5Meth_count.txt intron_3Meth_count.txt intron_5Meth_count.txt|awk '$2==0{print $0}' |sort |uniq -c|awk '$1==4{print $2"\t"0}' >No_Methlation_exon.txt
for i in 1
do
cat No_Methlation_exon.txt exon_3Meth_count.txt |sort |uniq -u >exon_3Meth_No0count.txt
cat No_Methlation_exon.txt exon_5Meth_count.txt |sort |uniq -u >exon_5Meth_No0count.txt
cat No_Methlation_exon.txt intron_3Meth_count.txt |sort |uniq -u >intron_3Meth_No0count.txt
cat No_Methlation_exon.txt intron_5Meth_count.txt |sort |uniq -u >intron_5Meth_No0count.txt
done
同样没有发现想要的效果,比如Alter 的剪切位点中exon和intron之间的甲基化程度相差不大,而Constitutive 剪切位点附近的甲基化程度相差很大

不把intron与exon拆开看了,就分析5'|3’剪切位点附近甲基化程度和CG含量
## 将同一个外显子附近的甲基化数据合并
cat *3Meth_No0count.txt|awk '{a[$1]+=$2}END{for(i in a) print "ConExon\t""3_splice\t"a[i]}' >3Meth_No0count.txt
cat *5Meth_No0count.txt|awk '{a[$1]+=$2}END{for(i in a) print "ConExon\t""5_splice\t"a[i]}' >5Meth_No0count.txt
佛了,想要差异大一些又没差异

参考
bedtools使用 https://www.jianshu.com/p/6c3b87301491
ggplot2 https://blog.csdn.net/Bone_ACE/article/details/47427453# 标签修改
Last updated
Was this helpful?