PNAS Fig.1的复现（续）

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)]

友情宣传

• 全国巡讲全球听（买一得五），第二期 ，你的生物信息学入门课
• 生信技能树的2019年终总结 ，你的生物信息学成长宝藏
• 2020学习主旋律，B站74小时免费教学视频为你领路