ChIP-Seq文献数据重新分析解读第一例

文章是:Genome-wide maps of H3K4me2/3 in prostate cancer cell line LNCaP,数据在GEO可以下载。GSE20042,下面的所有分析,需要26G的空间。
作者想看看用 dihydrotestosterone (雄激素)处理了 cancer cell line LNCaP 这个细胞系之后,看看组蛋白甲基化修饰变化,主要是看H3K4me2和H3K4me3这两种组蛋白甲基化区别,分成三组,分别是处理前,处理后4H和16H,共有5个条件的数据,但是有7个fastq文件。
测序仪是:Illumina Genome Analyzer (Homo sapiens)
主要是为了分析差异 核小体定位点 区别:Model for identifying differential transcription factor binding locations
作者在这里进行数据分析软件(NPS)很旧了,也是哈佛刘小乐实验室出品的,我这里就不用了。
数据处理详情如下:
Bed: Sequence reads were obtained and mapped to the human genome (March, 2006) using the Illumina Genome Analyzer Pipeline.
Peaks: Peak detection was performed with the "Nucleosome Positioning from Sequencing (NPS)" algorithm (http://liulab.dfci.harvard.edu/NPS/)
Processed data file build: hg18
所以我重新重复这 个数据分析,用的hg19,还有MACS2z这个软件
作者同时也测了芯片数据:Affymetrix U133 Plus 2.0 microarray data,但是似乎并没有给地址,我们先不管
首先下载数据
cd ~/CHIPseq_test/
mkdir GSE20042_H3K4me2_3 && cd GSE20042_H3K4me2_3
mkdir rawData && cd rawData
GSM503903    H3K4me2_Vehicle_ChIPSeq    SRR037146(Read 7,530,267 spots )/SRR037147(Read 6,215,981 spots )
GSM503904    H3K4me2_DHT_4h_ChIPSeq    SRR037148(Read 6,510,159 spots )/SRR037149(Read 6,246,716 spots )
GSM503905    H3K4me2_DHT_16h_ChIPSeq    SRR037150    Read 9,685,845 spots
GSM503906    H3K4me3_Vehicle_ChIPSeq    SRR037151    Read 6,755,854 spots
GSM503907    H3K4me3_DHT_4h_ChIPSeq    SRR037152    Read 4,761,769 spots
## 可以看到测序量并不大,因为文章比较老了,其实现在一般要测20M的reads
ls *sra |while read id; do ~/biosoft/sratoolkit/sratoolkit.2.6.3-centos_linux64/bin/fastq-dump $id;done
rm *sra
ls *.fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;done
mkdir QC_results
mv *zip *html QC_results/
##接下来做比对
## cat >run_bowtie2.sh  运行这个脚本批量做alignment
ls *.fastq | while read id ;
do
echo $id
~/biosoft/bowtie/bowtie2-2.2.9/bowtie2 -3 5 -p 8 -x ~/biosoft/bowtie/hg19_index/hg19 -U $id   -S ${id%%.*}.sam  2>${id%%.*}.align.log;
samtools view -bhS -q 30  ${id%%.*}.sam > ${id%%.*}.bam  ## -F 1548 https://broadinstitute.github.io/picard/explain-flags.html
samtools sort   ${id%%.*}.bam ${id%%.*}.sorted  ## prefix for the output
samtools index ${id%%.*}.sorted.bam
done
然后下载GEO的核小体定位点(peaks)结果:
tar xvf GSE20042_RAW.tar
ls *gz |xargs gunzip
wc -l *txt
  235639 GSM503903_LNCaP_H3K4me2_Vehicle_2Lanes_normalized_peak.txt
  248570 GSM503904_LNCaP_H3K4me2_DHT_4h_2Lanes_normalized_peak.txt
  185892 GSM503905_LNCaP_H3K4me2_DHT_16h_peak.txt
   74491 GSM503906_LNCaP_H3K4me3_Vehicle_normalized_peak.txt
  104022 GSM503907_LNCaP_H3K4me3_DHT_4h_normalized_peak.txt
然后根据比对的bam文件来可视化这些核小体peaks,很诡异,不知道他是如何找到的这些peaks,这些peaks画图之后根本看不出来,后来我才知道,是因为peaks的位点是hg18的坐标,而我用的是自己的bam文件来画图,所以~~~~
画图代码如下:
Rscript ~/CHIPseq_test/peakView.R GSM503907_LNCaP_H3K4me3_DHT_4h_normalized_peak.bed  ../rawData/SRR037152.sorted.bam
这个peakView.R代码很简单,就是用samtools depth命令提取每个peaks区域的坐标,然后画曲线即可
然后我用MACS2软件来call peaks 看看:
ls *sorted.bam |while read id;do ( nohup time ~/.local/bin/macs2 callpeak -t $id -f BAM -g hs -n ${id%%.*} 2>${id%%.*}.masc2.log &) ;done  
## 这里批量对7个测序文件做peaks callling 
mkdir ../MACS2results
mv *bed *xls *Peak *r  ../MACS2results
cd ../MACS2results
ls *.xls | while read id ;
do
echo $id
grep '^chr\S' $id |perl -alne '{print "$F[0]\t$F[1]\t$F[2]\t$F[9]\t$F[7]\t+"}' >${id%%.*}.bed
done
然后重新浏览peaks
Rscript ~/CHIPseq_test/peakView.R SRR037152_peaks.bed  ../rawData/SRR037152.sorted.bam
看起来我call的peaks还挺靠谱的, 图片以后再上传!

 

Comments are closed.