github：https://github.com/scottzijiezhang/m6A-seq_analysis_workflow
网页：https://scottzijiezhang.github.io/m6A-seq_analysis_workflow/
R包：scottzijiezhang/m6Amonster: Analyze m6A seq data version 0.1.2 from GitHub
R包的github：https://github.com/scottzijiezhang/m6Amonster
比对成bam
install.packages("devtools")
library(devtools)
install_github("scottzijiezhang/m6Amonster")
装个包真难！
ref:peakDistribution:
https://rdrr.io/github/scottzijiezhang/m6Amonster/src/R/peakDistribution.R
library(m6Amonster)
samplenames.a = c("DMSO1","DMSO2")
samplenames.b = c("OHT1","OHT2")
monster.a <- countReads(samplenames = samplenames.a,
gtf = "mm10.gtf",
bamFolder = "bam_files_dmso",
outputDir = "m6Amonster",
modification = "m6A",
binSize = 50,
threads = 20, paired=TRUE )
monster.b <- countReads(samplenames = samplenames.b,
gtf = "mm10.gtf",
bamFolder = "bam_files_oht",
outputDir = "m6Amonster",
modification = "m6A",
binSize = 50,
threads = 20 , paired=TRUE )
head(monster.a$reads)
head(monster.b$reads)
# 1.这里报错：chr的问题，把gtf里的去掉即可。
# 2.bam可以提前index；about 1 hour
monsterpeaks.a <- callPeakFisher(monster.a, min_counts = 15, peak_cutoff_fdr = 0.05 , peak_cutoff_oddRatio = 1, threads = 20)
monsterpeaks.b <- callPeakFisher(monster.b, min_counts = 15, peak_cutoff_fdr = 0.05 , peak_cutoff_oddRatio = 1, threads = 20)
【1 hour】
head(monsterpeaks.a$peakCallResult)
head(monsterpeaks.b$peakCallResult)
Joint_peak.a <- reportJointPeak(monsterpeaks.a, threads = 20)
Joint_peak.b <- reportJointPeak(monsterpeaks.b, threads = 20)
# uniq-peaks：
unique.Joint_peak.a <- Joint_peak.a[which(!duplicated(paste(Joint_peak.a$chr,Joint_peak.a$start,Joint_peak.a$end,sep = ":"))),]
unique.Joint_peak.b <- Joint_peak.b[which(!duplicated(paste(Joint_peak.b$chr,Joint_peak.b$start,Joint_peak.b$end,sep = ":"))),]
# 输出：
write.table(unique.Joint_peak.a, file = "dmso_Joint_peak.bed", sep = "\t", row.names = F, col.names = F, quote = F)
write.table(unique.Joint_peak.b, file = "oht_Joint_peak.bed", sep = "\t", row.names = F, col.names = F, quote = F)
# 画图：
plotMetaGene(Joint_peak.a,gtf = "mm10.gtf")
plotMetaGene(Joint_peak.b,gtf = "mm10.gtf")
# motif：
fa=/media/shen/6a524d78-97d1-481c-b068-8116a4d007f8/sun/refdata/gencode_GRCm38/GRCm38.chr.fa
bedtools getfasta -fi $fa -bed dmso_Joint_peak.bed -fo dmso_Joint_peak.fa -split -s
bedtools getfasta -fi $fa -bed oht_Joint_peak.bed -fo oht_Joint_peak.fa -split -s
findMotifs.pl dmso_Joint_peak.fa fasta homer_dmso_motif -fasta ran.fa -rna -p 20 -len 5,6,7
findMotifs.pl oht_Joint_peak.fa fasta homer_oht_motif -fasta ran.fa -rna -p 20 -len 5,6,7
下一步分析安装R包：RADAR：
install.packages("devtools")
library(devtools)
install_github("scottzijiezhang/RADAR")
library("RADAR")
