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%

  1. 四倍体陆地棉

    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%

  2. 二倍体雷蒙德氏棉

    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# 标签修改

Previous20200102GO富集分析Next甲基化重复间的处理

Last updated