methy pipeline
2024.11.26更新
https://github.com/jsh58/DMRfinder
- 将不管从什么方法得到的甲基化文件转换成如下格式：
qury.methy.bed: （pbmm2等方法得到的甲基化位点的bed文件）
chr1 3050095 3050096 100 8 0
chr1 3050096 3050097 83.3 5 1
...
具体含义为：
染色体 起始位置 起始位置+1 甲基化水平（100分制） 甲基化的C的数量 未甲基化的C的数量
- 将qury甲基化区域转换为参考（mouse）的坐标
注，在这里，已经将小鼠的甲基化区域处理好并转换成以下格式：[mm_methy.region.bed]
mm_chr1 3050095 3050404 mm_chr1_3050095_3050404
mm_chr1 3052296 3052381 mm_chr1_3052296_3052381
mm_chr1 3052735 3052828 mm_chr1_3052735_3052828
...
具体含义为 染色体 起始位置 终止位置 区域标记
liftOver mm_methy.region.bed target.qury.all.chain qury2target_methy.map.loc qury2target_methy.unmap.loc \
-minMatch=0.25
- 获取qury的甲基化区域
bedtools sort -i qury2target_methy.map.loc > qury2target_methy.map.loc.sort.bed
bedtools intersect -a qury2target_methy.map.loc.sort.bed -b qury.methy.clean.bed -wo > test.inter
其中：qury.methy.clean.bed（与step1中含义相同）
carmeli_Chr1 2902 2903 1.0 1 0
carmeli_Chr1 2924 2925 1.0 1 0
...
染色体 起始位置 起始位置+1 甲基化水平（1） 甲基化的位点 未甲基化的位点
python3 01.py test.inter > qury2target.conbineCpG.bed
bedtools sort -i qury2target.conbineCpG.bed > qury2target.conbineCpG.sort.bed
python3 02.py target.conbineCpG.bed qury2target.conbineCpG.sort.bed > 01target_qury.conbineCpG.input
根据https://github.com/jsh58/DMRfinder中添加列名：例如，
chr start end CpG carmeli-N carmeli-X mm-N mm-X
mm_chr1 3284715 3284827 4 60 58 60 42
mm_chr1 3470125 3470230 4 60 52 71 18
...
python3 03.py 01target_qury.conbineCpG.input > 01target_qury.conbineCpG.filter.input
Rscript ~/software/DMRfinder/findDMRs.v4.r -i 01mm_carmeli.conbineCpG.final.filter.input -o test.output \
-n Sca,mouse carmeli mm
galili pipeline
ccsmethy结果处理
cat galili.methy.bed| cut -f 1-3,10,11 | awk '{print$1"\t"$2"\t"$3"\t"$4"\t"$4*$5/100"\t"$4-$4*$5/100"\t"$5}' > ../01.cleandata/galil
i.methy.bed
1. 提取全基因组各区域位置
1.1 处理重复序列文件
只提取TE的bed文件
python3 00.filter.repeat.py galili.repeat.gff > galili.TE.bed
提取重复序列位置的bed文件
cat galili.repeat.gff | cut -f 1,4,5,9 | awk '{print$1"\t"$2-1"\t"$3"\t"$4}' | bedtools sort -i - > galili.repeat.sort.bed
1.2 处理基因区位置文件
基因区
cat galili.gene.gff | awk '$3=="gene"{print$0}' | cut -f 1,4,5 | awk '{print$1"\t"$2-1"\t"$3}' > galili.gene.bed
外显子区
cat galili.gene.gff | awk '$3=="exon"{print$0}' | cut -f 1,4,5 | awk '{print$1"\t"$2-1"\t"$3}' > galili.exon.bed
内含子区
bedtools subtract -a galili.gene.bed -b galili.exon.bed > galili.intron.bed
1.3 基因间区
cat galili.gene.bed galili.repeat.sort.bed | bedtools sort -i - > galili.function.sort.bed
cat Galili.Chr.v3.fasta.fai | awk '{print$1"\t"$2"\t"$2}' | bedtools sort -i - | cut -f 1-2 > galili.genome.len
bedtools complement -i galili.function.sort.bed -g galili.genome.len > galili.intergentic.bed
1.4 计算各区域甲基化含量
bedtools intersect -a galili.exon.bed -b galili.methy.bed -wb > galili.exon.methy.output
bedtools intersect -a galili.intergentic.bed -b galili.methy.bed -wb > galili.intergentic.methy.output
bedtools intersect -a galili.intron.bed -b galili.methy.bed -wb > galili.intron.methy.output
bedtools intersect -a galili.repeat.sort.bed -b galili.methy.bed -wb > galili.repeat.methy.output
bedtools intersect -a galili.TE.bed -b galili.methy.bed -wb > galili.TE.methy.output
1.5 转换成画图格式
cat galili.exon.methy.output | cut -f 10 | awk '{print$1"\texon\tgalili"}' > galili.exon.methy.output.plot
cat galili.intergentic.methy.output | cut -f 10 | awk '{print$1"\tintergentic\tgalili"}' > galili.intergentic.methy.output.plot
cat galili.intron.methy.output | cut -f 10 | awk '{print$1"\tintron\tgalili"}' > galili.intron.methy.output.plot
cat galili.repeat.methy.output | cut -f 10 | awk '{print$1"\trepeat\tgalili"}' > galili.repeat.methy.output.plot
cat galili.TE.methy.output | cut -f 11 | awk '{print$1"\tTE\tgalili"}' > galili.TE.methy.output.plot
cat *.plot > galili.genomic.in.plot
around TSS
1. 提取转录起始位点
cat ../00.rawdata/galili.gene.gff | grep -w "exon1" | awk '{print$1"\t"$4"\t"$4+1}' > galili.TSS.bed
2. 提取甲基化位点到TSS的距离
bedtools closest -a ../../01.cleandata/galili.methy.bed -b galili.TSS.bed -D "b" > dist_TSS.bed
3. 过滤上下游3kb
cat dist_TSS.bed | cut -f 7,11 | awk '$2<3000{print$0}' | awk '$2>-3000{print$0}' > dist_TSS_in3K.plot
4. 画图（R）
library(tidyverse)
library(ggplot2)
library(dplyr)
data <- read.table("dist_TSS_in3K.plot",header = F,
col.names = c("me_level","distance")) %>%
mutate(other="kda") %>%
select(other,me_level,distance)
data <- data[which(data$distance < 3000 & data$distance > -3000),2:3]
data <- data[order(data[,2]),]
results_dist <- apply(data,2,function(x) ave(x, data[,2], FUN=mean))
results_dist <- results_dist[!(duplicated(results_dist[,2])),]
plot(x=results_dist[,2], y=results_dist[,1], type="n", ylim=c(0,100),
xlab="Dist to TSS", ylab="Average Methylation")
fit <- lm(results_dist[,1]~poly(results_dist[,2],6,raw=T))
lines(results_dist[,2],predict(fit,data.frame(x=results_dist[,1])),
lwd=3, col="black")
around TES
1. 提取转录终止位点
cat ../..//00.rawdata/galili.gene.gff | awk '$3=="exon"{print$0}' > galili.exon.locate
python3 01.extract.TES.py galili.exon.locate > galili.TES.bed
2. 提取甲基化位点到TES的距离
bedtools closest -a ../../01.cleandata/galili.methy.bed -b galili.TES.bed -D "b" > dist_TES.bed
3. 过滤上下游3kb
cat dist_TES.bed | cut -f 7,11 | awk '$2<3000{print$0}' | awk '$2>-3000{print$0}' > dist_TES_in3K.plot
4. 画图
基因体区甲基化差异
cat scg.txt | cut -f 4 | tr "," "\t" > species_ortho.txt
cat species_ortho.txt | cut -f 1 > mouse.gene.list
cat species_ortho.txt | cut -f 2 | sed 's/model/TU/'> galili.gene.list
cat galili.final.gff | awk '$3=="gene"{print$0}' | cut -f 1,4,5,9 | tr ";" "\t" | cut -f 1-4 | sed 's/ID=//' > galili.gene.locat.txt
cat mus.final.gff | awk '$3=="mRNA"{print$0}' | cut -f 1,4,5,9 | tr ";" "\t" | cut -f 1-4 | sed 's/ID=//' > mouse.gene.locat.txt
python3 00.py galili.gene.locat.txt galili.gene.list > galili.ortho.gene.locate.bed
bedtools sort -i galili.ortho.gene.locate.bed > galili.ortho.gene.locate.sort.bed
python3 02.py mouse.ortho.gene.locate.sort.bed mouse.methy.new.bed 01.py mouse > mouse.run.sh
python3 02.py galili.ortho.gene.locate.sort.bed galili.methy.new.bed 01.py galili > galili.run.sh
sh mouse.run.sh && sh galili.run.sh
cat galili/*.stat > galili.ortho.gene.methy_level.output
cat mouse/*.stat > mouse.ortho.gene.methy_level.output
python3 03.py species_ortho.rename.txt mouse.ortho.gene.methy_level.output galili.ortho.gene.methy_level.output > species_ortho.methy.compart.stat
cat species_ortho.methy.compart.stat | grep -v "p_value" | awk '$1<=0.05{print$0}' | awk '$2>=2{print$0}' | cut -f 3 | sed 's/rna-//' > galili.hypomethy.mRNA_symbol.list
cat mus.final.gff | awk '$3=="mRNA"{print$0}' | awk '{print$9}' | tr ";" "\t" | cut -f 1,2 | sed 's/ID=rna-//' | sed 's/Parent=gene-//' > mouse.gene.duizhao.list
python3 04.py mouse.gene.duizhao.list galili.hypomethy.mRNA_symbol.list > galili.hypomethy.gene_symbol.list
基因启动子区甲基化差异
cat mus.final.gff | awk '$3=="exon"{print$0}' | grep "\-1\;" | awk -F "\t" '{print$1"\t"$4"\t"$5"\t"$9}' | tr ";" "\t" | cut -f 1-3,5 | sed 's/Parent=//' | awk '{print$1"\t"$2-2000"\t"$2+50"\t"$4}' > mouse.prometer.locate.txt
cat galili.final.gff | grep -w "exon1" | awk -F "\t" '{print$1"\t"$4"\t"$5"\t"$9}' | tr ";" "\t" | cut -f 1-3,5 | sed 's/Parent=//' | awk '{print$1"\t"$2-2000"\t"$2+50"\t"$4}' > galili.prometer.locate.txt
【用sed把-值改为0】
python3 00.py galili.prometer.locate.txt galili.gene.list > galili.ortho.promoter.locate.bed
bedtools sort -i galili.ortho.promoter.locate.bed > galili.ortho.promoter.locate.sort.bed
python3 00.py mouse.prometer.locate.txt mouse.gene.list > mouse.ortho.promoter.locate.bed
bedtools sort -i mouse.ortho.promoter.locate.bed > mouse.ortho.promoter.locate.sort.bed
python3 02.py galili.ortho.promoter.locate.sort.bed galili.methy.bed 01.py galili > galili.run.sh
验证TSS双峰
cat galili.prometer.locate.txt | awk '{print$1"\t"$2-1"\t"$3"\t"$4}' > galili.prometer.locate.bed
cat mouse.prometer.locate.txt | awk '{print$1"\t"$2-1"\t"$3"\t"$4}' > mouse.prometer.locate.bed
sed -i 's/-[0-9]\+/0/' mouse.prometer.locate.bed
sed -i 's/-[0-9]\+/0/' galili.prometer.locate.bed
bedtools intersect -a galili.prometer.locate.bed -b galili.methy.bed -wb > galili.promoter.methyLevel.output
bedtools intersect -a mouse.prometer.locate.bed -b mouse.methy.new.bed -wb > mouse.promoter.methyLevel.output
cat galili.promoter.methyLevel.output | awk '{print"galili\t"$11}' > galili.promoter.plot
cat mouse.promoter.methyLevel.output | awk '{print"mouse\t"$8}' > mouse.promoter.plot
