【基本流程】Bismark—Methylkit—Genomation—GO/KEGG
只求自己记得 如果有人搜到了这个 请你只看上面的四个基本流程 然后分别搜索他们的用法
比较坑的地方在于用Genomation包+BED文件来注释,以及将ID转换为NAME进而用于GO/KEGG
https://blog.csdn.net/qq_29300341/article/details/80551131
(btw:这篇文章写得不错)
bismark将bismark的结果导入R,用R包methylKit进行下一步分析
这里需要将结果转换一下格式
使用的文件是*.CpG_report.txt(也就是下文加粗的那个)自己用时自己替换
awk ' { if ( $4+$5!=0 && $3~/+/) print $1"."$2,$1,$2,"F",$4+$5,sprintf("%.2f", $4/($4+$5)),sprintf("%.2f", 100* $5/($4+$5)); else if ($4+$5!=0 && $3~/-/) print $1"."$2,$1,$2,"R",$4+$5,sprintf("%.2f", 100*$4/($4+$5)),sprintf("%.2f", $5/($4+$5)) }' OFS="\t" M4-WYWPE19011596_R1_val_1_bismark_bt2_pe.CpG_report.txt > tmp.txt && sed "1ichrBase\tchr\tbase\tstrand\tcoverage\tfreqC\tfreqT" tmp.txt > M4.CpG.pro.txt && rm tmp.txt
另一种方法
然后就可以导入R了
接着是下载BED文件用来注释
UCSC的Table Browser
http://genome-asia.ucsc.edu/cgi-bin/hgTables?hgsid=730113983_azih0C4PgtTEjwhkyyid7GUaRnj0&clade=mammal&org=Mouse&db=mm10&hgta_group=genes&hgta_track=refSeqComposite&hgta_table=0&hgta_regionType=genome&position=chr12%3A56%2C694%2C976-56%2C714%2C605&hgta_outputType=primaryTable&hgta_outFileName=name.txt
然后用genomation注释
将Gene ID转换为Gene Symbol
ENTREZID<- bitr(a[,4], fromType = "ENSEMBLTRANS", toType="SYMBOL",OrgDb = org.Mm.eg.db)
把Gene Symbol用于GO/KEGG
网友评论