1. 写在前面

参考的技能树视频链接：https://www.bilibili.com/video/BV1KJ411X7Zx?t=2&p=22
视频用到的文章的部分图表是我之前复现过的，笔记链接：https://www.jianshu.com/p/d78121da31c4
故本文仅补充一些上次复现时没有学到的小知识

2. 通过 bw 文件计算 ChIP-seq 样本间相关性

两个命令均来自 deeptools
嫌丑就载入R调整

multiBigwigSummary bins -b ../align/bla.*.bw -o results.npz -p 24
plotCorrelation -in results.npz \
--corMethod spearman --skipZeros \
--plotTitle "bla" \
--whatToPlot heatmap --colorMap RdylBu --plotNumbers \
--plotFileFormat pdf \
-o bla.pdf \
-- outFileCorMatrix SpearmanCorr_readCounts.tab

3. 批量化操作的一个思路

# 别整天for循环了
anno_bed <- function(bedPeaksFile){}
anno_result = lapply(list.files(path = '', pattern = '', full.names = T), anno_bed)

4. 转换 bed 文件的基因组版本

以 dm3 转 dm6 为例
方法一：使用 UCSC 的 liftover 网页工具：https://genome.ucsc.edu/cgi-bin/hgLiftOver

方法二：使用 R 包rtracklayer,liftOver

首先需要一个坐标注释文件：https://hgdownload-test.cse.ucsc.edu/goldenPath/dm3/liftOver/，名称为 dm3ToDm6.over.chain.gz
代码

library(rtracklayer)
library(liftOver)
path = 'bla/dm3ToDm6.over.chain'
ch = import.chain(path)
ch
library(ChIPseeker)
bed = readPeakFile('bla/bla.bed')
seqlevelsStyle(bed) = "UCSC"  # necessary
new_bed = liftOver(bed, ch)
class(new_bed)
# 还需要研究一下怎么把 GRanges 对象写出为 bed 文件

5. 获得差异表达基因的坐标

library(TxDb.Dmelanogaster.UCSC.dm6.ensGene)
txdb = TxDb.Dmelanogaster.UCSC.dm6.ensGene
gene_cor = as.data.frame(genes(txdb))
library(org.Dm.eg.db)
fb_table = toTable(org.Dm.egFLYBASE)
symbol_table = toTable(org.Dm.egSYMBOL)
fb_to_symbol = merge(fb_table,symbol_table,by = 'gene_id')
colnames(gene_cor)[6] = 'flybase_id'
gene_cor_symbol = merge(gene_cor,fb_to_symbol,by = 'flybase_id')
goi_cor = gene_cor_symbol[match(goi[goi %in% gene_cor_symbol$symbol],
                                gene_cor_symbol$symbol),c(4:6)]