前言

最近一直在看小麦表观遗传相关的文献，然后下载了文章中的数据进行分析，分析的流程大致跟网上教程相似，但是由于小麦基因组比较大（16G），研究不像人类，小鼠及一些模式植物广泛，完善，所以中间有些分析内容还是有些不同，于是把前期分析的代码放在这里，以供交流学习
数据出处：https://doi.org/10.1186/s13059-020-01998-1

质控比对

按照文章中的标准进行过滤和比对，因为样品比较少，所以没有批处理，也更直观一些，中间可能挑了单个样品跑了测试，所以有的样品没有出现在主脚本中

文章处理流程

#QC trim
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K27me3_1.fastq.gz H3K27me3_2.fastq.gz H3K27me3_1_trimmed.fastq.gz H3K27me3_1_unpaired.fastq.gz H3K27me3_2_trimmed.fastq.gz H3K27me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K36me3_1.fastq.gz H3K36me3_2.fastq.gz H3K36me3_1_trimmed.fastq.gz H3K36me3_1_unpaired.fastq.gz H3K36me3_2_trimmed.fastq.gz H3K36me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K4me3_1.fastq.gz H3K4me3_2.fastq.gz H3K4me3_1_trimmed.fastq.gz H3K4me3_1_unpaired.fastq.gz H3K4me3_2_trimmed.fastq.gz H3K4me3_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 H3K9ac_1.fastq.gz H3K9ac_2.fastq.gz H3K9ac_1_trimmed.fastq.gz H3K9ac_1_unpaired.fastq.gz H3K9ac_2_trimmed.fastq.gz H3K9ac_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
java -jar /usrdata/users/hwwang/software/Trimmomatic-0.36/trimmomatic-0.36.jar PE -phred33 -threads 16 RNAPII_1.fastq.gz RNAPII_2.fastq.gz RNAPII_1_trimmed.fastq.gz RNAPII_1_unpaired.fastq.gz RNAPII_2_trimmed.fastq.gz RNAPII_2_unpaired.fastq.gz -validatePairs ILLUMINACLIP:/usrdata/users/hwwang/software/Trimmomatic-0.36/adapters/TruSeq3-PE.fa:2:30:10 LEADING:5 TRAILING:5 MINLEN:20 2>>trim.log
#align
ref=/usrdata/users/hwwang/hychao/genome/iwgsc_refseqv1.0_all_chromosomes/bowtie2_index/wheat_bowtie_index
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K27me3_1_trimmed.fastq.gz -2 H3K27me3_2_trimmed.fastq.gz > H3K27me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K36me3_1_trimmed.fastq.gz -2 H3K36me3_2_trimmed.fastq.gz > H3K36me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K4me3_1_trimmed.fastq.gz -2 H3K4me3_2_trimmed.fastq.gz > H3K4me3.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 H3K9ac_1_trimmed.fastq.gz -2 H3K9ac_2_trimmed.fastq.gz > H3K9ac.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 Input_1_trimmed.fastq.gz -2 Input_2_trimmed.fastq.gz > Input.sam
bowtie2 -p 12 --very-sensitive -x $ref -1 RNAPII_1_trimmed.fastq.gz -2 RNAPII_2_trimmed.fastq.gz > RNAPII.sam
#sort
samtools sort -@ 12 -o H3K27me3.bam H3K27me3.sam
samtools sort -@ 12 -o H3K36me3.bam H3K36me3.sam
samtools sort -@ 12 -o H3K4me3.bam H3K4me3.sam
samtools sort -@ 12 -o H3K9ac.bam H3K9ac.sam
samtools sort -@ 12 -o Input.bam Input.sam
samtools sort -@ 12 -o RNAPII.bam RNAPII.sam
#filter
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K27me3.bam > H3K27me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K36me3.bam > H3K36me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K4me3.bam > H3K4me3_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' H3K9ac.bam > H3K9ac_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' Input.bam > Input_aln.bam
sambamba view -h -t 8 -f bam -F 'not unmapped and not duplicate and mapping_quality >= 10' RNAPII.bam > RNAPII_aln.bam

生成比对结果

下面需要注意的一点就是samtools建立index生成 .bai 文件对染色体的长度是有限制的，因为小麦染色体长度在7，800M左右，运行的时候会报错，所以在后面加一个-c参数，生成.csi文件

#qc_bam
ls  *aln.bam  |xargs -i samtools index -c {} 
ls  *aln.bam  | while read id ;do (nohup samtools flagstat $id > $(basename $id "aln.bam").stat & );done
grep N/A *.stat|grep %

IGV可视化

下一步就是生成IGV可视化的文件，这一步有一点小坑，就是bamCoverage软件在进行归一化的时候有一点改动，下面是之前的代码

#igv calculate reads density
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K27me3_aln.bam -o H3K27me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K36me3_aln.bam -o H3K36me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K4me3_aln.bam -o H3K4me3.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b H3K9ac_aln.bam -o H3K9ac.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b Input_aln.bam -o Input.bigwig
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeTo1x 14547261565 -b RNAPII_aln.bam -o RNAPII.bigwig

这里是我修改的

bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K36me3_aln.bam -o H3K36me3.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K4me3_aln.bam -o H3K4me3.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b H3K9ac_aln.bam -o H3K9ac.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b Input_aln.bam -o Input.bigwig 2>>run.log
bamCoverage -p 8 -bs 100 -of bigwig -extendReads --normalizeUsing RPGC --effectiveGenomeSize 14547261565 -b RNAPII_aln.bam -o RNAPII.bigwig 2>>run.log

总的来说，跟文章中的处理过程是一致的

文章处理流程

最后就是call peak了！

#call peaks(histone)
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K27me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K36me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K4me3_result --outdir  ./
macs2 callpeak -f BAM --nomodel --to-large -p 0.01 -broad -g 17e9 --bw 300 -c Input_aln.bam -t -n H3K9ac_result --outdir  ./
#call peaks RNAPII
macs2 callpeak -f BAMPE --nomodel –q 0.001 --broad-cutoff 0.01 -g 17e9 --bw 300 -c Input_aln.bam -t -n RNAPIIRNAPII_result --outdir ./

上面的call peak代码根据文章所述

文章分析流程

后面的分析内容先简单放点吧，因为还没有全部完成，挑了一个数据，做了peak注释

setwd("E:\\R_file\\test_chip\\")
library("ChIPseeker")
library("GenomicFeatures")
wheat_txdb <- loadDb("E:\\R_file\\wheat.sqlite")
RNAPII <- readPeakFile("H3K27me3_result_peaks.broadPeak.fold3")
peakAnno <- annotatePeak(RNAPII,tssRegion=c(-3000, 3000),TxDb=wheat_txdb)
plotAnnoPie(peakAnno)
plotAnnoBar(peakAnno)
vennpie(peakAnno)
upsetplot(peakAnno)

随便几个图

更新一下

输出注释的peaks，一些图可以尝试画一下

png("covplot_H3K27me3.png")
covplot(H3K27me3,weightCol=5)
dev.off()
##get promoter region +-3Kb and write out
promoter <- getPromoters(TxDb=txdb, upstream=3000, downstream=3000)
write.table(promoter,"wheat_promoter.txt",row.names=FALSE,sep="\t")
#tagMatrix <- getTagMatrix(H3K27me3, windows=promoter)
#png("tagHeatmap.png")
#tagHeatmap(tagMatrix, xlim=c(-3000, 3000), color="red")
#dev.off()

##peak anno
peakAnno <- annotatePeak(H3K27me3,tssRegion=c(-3000, 3000),TxDb=txdb)
anno_out=as.data.frame(peakAnno)
write.table(anno_out,"H3K27me3_peakAnno.txt",row.names=FALSE,sep="\t")

标记数目比较多的话，可以使用 chromHMM 软件分析不同标记组合的染色质状态

10种states

distribution

写在最后

总的来说，前面的分析并不麻烦，就是有一些地方容易踩坑，写下这些，主要是为了方便回看学习，也为了后面有做这方面的同学不要再犯跟我一样的错误。关于小麦的 txDb 文件，怎样制作网上都有教程，搜一下就可以了