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小时免费教学视频为你领路
网友评论