很久不见。
错过了很多的询问和回复,很抱歉。
这段时间在整理我的第一篇paper,大部分涉及的东西不太具有可分享性(其实也是应该记录的)
———————————————————————————————
今天的笔记是记录处理ENCODE来源的ChIP-seq下载及分析的过程。
———————————————————————————————
1. ENCODE是什么?
用了别人的数据一定要表示尊重的。
ENCODE的全称是The Encyclopedia of DNA Elements (DNA元件百科全书),是一个由美国人类基因组研究所在2003年9月发起至今仍在维护和更新的一个公共联合研究项目(https://www.encodeproject.org/)。
ENCODE现有的数据主要分为了三大块:
(1) Functional genomics: 这部分的数据是我们最常会用到的数据,包括的实验来源数据有
① TF ChIP-seq; Histone ChIP-seq; DNase-seq; Mint-ChIP-seq; FAIRE-seq; MNase-seq; ATAC-seq; snATAC-seq; DNase-HS; DNAme array; WGBS; RRBS; MeDIP-seq; Hi-C; intact Hi-C; in situ Hi-C; ChIA-PET; Repli-seq; Repli-chip; PAS-seq; WGS;
② total RNA-seq; polyA plus RNA-seq; microRNA-seq; scRNA-seq; small RNA-seq; long read RNA-seq; RAMPAGE; RNA microarray; genotyping array; CAGE; microRNA counts; 5‘RLM RACE;
③ eCLIP; iCLIP; RNA Bind-n-seq; RIP-chip; RIP-seq; PRO-cap; GRO-cap; PRO-seq; Circulome-seq等等...
(2) Functional characterization: 包括的数据有
enhancer reporter assay; Flow-FISH CRISPR screen; MPRA; proliferation CRISPR screen; STARR-seq; FACS CRISPR screen; perturbation followed by snATAC-seq; perturbation followed by scRNA-seq; CRISPR screen。
(3) Encyclopedia of elements: 这部分的数据经常联合(1)一起进行分析,包括imputation; candidate Cis-Regulatory Elements; chromatin state; representative DNase hypersensitivity sites。
除了这三个之外,ENCODE还有很多专项的研究计划,涵盖了人类样本human donor、类器官、细胞系等。详见主页。
ENCODE计划产出的数据为认识基因组功能元件提供了巨大的数据瑰宝。
2. ENCODE数据的下载
这里以下载Hela细胞ChIP-seq数据为例说明。
首先根据前面的介绍我们知道了ChIP-seq的数据都在Functional genomic这一部分里,点击进入这一部分之后,得到下面界面:
得到下载的数据可以通过左边的搜索框确定,也可以通过右边的matrix点击相应的格子进入。
这里我们点击了Hela S3对应的TF ChIP-seq方格之后进入以下界面:
在这里,左边的Assay, Analysis等信息的确定是我们下载的关键
对应的实验靶标,样本等根据自己的需求选择就行,这里需要说明的是一定不要忘记设置Analysis;
如图,Analysis里面的Available file types包括了fastq; bam; bed idr_ranked_peak; bed narrowPeak; bigBed narrowPeak; bigwig; bed broadPeak; bigBed broadPeak。真的很丰富!
在bam里,又有经过了picard等处理的alignment_bam以及unfiltered_bam。 也是根据自己的分析需求进行下载,不用重复下载耽误时间,详细文件说明参见https://www.encodeproject.org/chip-seq/transcription_factor/
都设定好之后,我们点击页面右边上面的Download,会出现以下信息:
这里,(1) Download default files是你经过筛选之后得到的想要下载文件的txt
(2)其中xargs 命令则是需要在Linux终端输入的根据下载文件进行数据下载所用到的语句。
xargs -L 1 curl -O -J -L < downloaded_files.txt
下载完成后的fastq、bam文件一定记得检查md5
3. bam文件的处理
在分析的过程中(1)需要进行特定区间的reads的提取; (2) 对bam进行filteration: 去除mismatch > 4的reads以及去除soft-clipped reads。
这部分是这个笔记的精华。也是找了很多资料之后柳岸花明终于实现的目前已知的对新手来说最能理解和操作的方法。
理解这两个需求我们需要复习bam文件
这部分已经有很多很多笔记和教程了,例如https://www.jianshu.com/p/7d15173540ae;大家自行观看。
要明确的是mismatch的信息储存在bam文件第12列的tag里,以XM表示。而soft-clipp reads则储存在bam的cigar里
—————————————————————————————
(1) 从bam文件中提取特定基因组区间的reads:
· 当我们需要提取一个区间时,可以使用smatools轻松达到目的:
samtools view -hb chr:start-end wgs.sort.bam > target.region.bam
· 当我们需要提取多个区间时,使用:
samtools view -hb -L target.bed deduped.bam > samtools_view_L_target.bam
更多的实现方式见:https://blog.csdn.net/tanzuozhev/article/details/88975801
(2) 去除mismatch > 4的reads;去除soft-clip reads:
这一部分,也找到了很多基于samtools + awk实现的方法,但是最终我使用的是nf-core/mnaseseq (https://hub.docker.com/r/nfcore/mnaseseq) 里面的一个脚本,命令如下:
bamtools filter -in test.bam -out test.filter.bam -script extraction.json
extraction.json的code如下:
{
"filters" : [
{ "id" : "mismatch",
"tag" : "NM:<4"
},
{ "id" : "cigar",
"cigar" : "*S*"
},
{ "id" : "insertion",
"cigar" : "*I*"
},
{ "id" : "deletion",
"cigar" : "*D*"
}
],
"rule" : " mismatch & !cigar & !insertion & !deletion "
}
bamtools 的使用说明参见:https://hcc.unl.edu/docs/applications/app_specific/bioinformatics_tools/data_manipulation_tools/bamtools/running_bamtools_commands/