文章是：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
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_3mkdir rawData && cd rawDatafor ((i=146;i<153;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP002/SRP002077/SRR037$i/SRR037$i.sra;done
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;donerm *srals *.fastq | while read id ; do ~/biosoft/fastqc/FastQC/fastqc $id;donemkdir QC_resultsmv *zip *html QC_results/##接下来做比对## cat >run_bowtie2.sh 运行这个脚本批量做alignmentls *.fastq | while read id ;doecho $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.htmlsamtools sort ${id%%.*}.bam ${id%%.*}.sorted ## prefix for the outputsamtools index ${id%%.*}.sorted.bamdone
然后下载GEO的核小体定位点(peaks)结果：
tar xvf GSE20042_RAW.tarls *gz |xargs gunzipwc -l *txt235639 GSM503903_LNCaP_H3K4me2_Vehicle_2Lanes_normalized_peak.txt248570 GSM503904_LNCaP_H3K4me2_DHT_4h_2Lanes_normalized_peak.txt185892 GSM503905_LNCaP_H3K4me2_DHT_16h_peak.txt74491 GSM503906_LNCaP_H3K4me3_Vehicle_normalized_peak.txt104022 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还挺靠谱的， 图片以后再上传！