搜索
查看: 4642|回复: 0

ChIP-Seq数据分析实战-H3K4me2和H3K4me3这两种组蛋白甲基化

[复制链接]

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2016-9-6 07:15:39 | 显示全部楼层 |阅读模式
作者想看看用 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,还有MACS2这个软件
作者同时也测了芯片数据:Affymetrix U133 Plus 2.0 microarray data,但是似乎并没有给地址,我们先不管
数据官网是:wide maps of H3K4me2/3 in prostate cancer cell line LNCaP http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE20042
首先也需要通读该文章,然后根据CHIP-seq数据分析流程来处理数据,验证作者的结果。

首先下载数据
cd ~/CHIPseq_test/
mkdir GSE20042_H3K4me2_3 && cd GSE20042_H3K4me2_3
mkdir rawData && cd rawData
for ((i=146;i<153;i++)) ;do wget ftp://ftp-trace.ncbi.nlm.nih.gov ... 02/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;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
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

然后我用MACS2软件来call peaks 看看:
# http://www2.uef.fi/documents/1698400/2466431/Macs2/f4d12870-34f9-43ef-bf0d-f5d087267602
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还挺靠谱的, 图片以后再上传!





上一篇:ChIP-Seq数据分析实战-喜树碱处理对TP53的结合能力的改变
下一篇:生信人必须了解的常见ftp资源中心
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-3-26 19:03 , Processed in 0.045163 second(s), 32 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.