全文分析流程学习按照：九月学徒ChIP-seq学习成果展
一、 怎么将SAR文件转为fastq文件？
1. 【方法】利用软件sratoolkit，下载后解压，找到其中的 fastq-dump.exe 可执行文件
（export PATH方法没用，用的alias方法。Export路径写到bin , alias路径写到可执行文件的绝对路径）
fastq-dump -I --split-3 SRR948800.sra -O /output/path/ #使用方法
#其中split-3表示如果是单端测序则一个sra文件出来一个fastq文件，如果是双末端，则一个sra问件对应两个fastq文件。
2. 【应用】写一个脚本：
Vim sra_to_fastq
#!bin/bash
# sra conver to fastq
#BSUB -J align145
#BSUB -n 10
#BSUB -R span[hosts=1]
#BSUB -o %J.out
#BSUB -e %J.err
#BSUB -q normal
nohup fastq-dump -I --split-3 /public/home/thu/chip_seq/sra/SRR*.sra &
bash sra_to_fastq
ps
##查看进程
Jobs -l
##查看nohup & 挂在后台的进程情况
二、 怎么用fastQC进行质控分析
【应用】写一个脚本：
#!bin/bash
# sra conver to fastq
#BSUB -J align145
#BSUB -n 10
#BSUB -R span[hosts=1]
#BSUB -o %J.out
#BSUB -e %J.err
#BSUB -q normal
nohup fastqc -o ./fastqc *.fastq &
三、 用multiqc整合报告
- multiqc *zip -o ./multiqc/
四、过滤低质量的fq文件
1. 【方法】用Trim Galore软件
说明书：Taking appropriate QC measures for RRBS-type or other -Seq applications with Trim Galore!
2. 【应用】运行一个脚本
1 #!bin/bash
2 #trim_galore
3 #BSUB -J align145
4 #BSUB -n 10
5 #BSUB -R span[hosts=1]
6 #BSUB -o %J.out
7 #BSUB -e %J.err
8 #BSUB -q normal
9
10 fq1=/public/home/thu/chip_seq/sra/SRR6795678_1.fastq.gz
11 fq2=/public/home/thu/chip_seq/sra/SRR6795678_2.fastq.gz
12
13 nohup trim_galore -q 20 --phred33 \
14 --paired \ #双端测序
15 --length 40 \ #小于40的read删除（可变，大小浮动对后续分析影响不大）
16 --fastqc \ #同时生成fastqc文件
17 --stringency 3 \
18 $fq1 $fq2 \ #fastq文件
19 --gzip & #压缩形式
3. trim_galore报告（若使用后台命令，则报告在nohup.out文件中）：
SUMMARISING RUN PARAMETERS
==========================
Input filename: ENCFF108UVC.fastq.gz #输入文件
Trimming mode: paired-end #双端测序结果
Trim Galore version: 0.5.0
Cutadapt version: 1.18
Quality Phred score cutoff: 25 #质量低于25的被删除，默认20（99%）
Quality encoding type selected: ASCII+33 #质量得分算法类型
Adapter sequence: 'AGATCGGAAGAGC' (Illumina TruSeq, Sanger iPCR; auto-detected)
Maximum trimming error rate: 0.1 (default) #最大错误率
Minimum required adapter overlap (stringency): 5 bp #最大adapter重复数目？
Minimum required sequence length for both reads before a sequence pair gets removed: 75 bp #最小的长度数量（可以适当放宽）
Output file will be GZIP compressed
This is cutadapt 1.18 with Python 2.7.6
Command line parameters: -f fastq -e 0.1 -q 25 -O 5 -a AGATCGGAAGAGC ENCFF108UVC.fastq.gz
Processing reads on 1 core in single-end mode ...
Finished in 1038.93 s (40 us/read; 1.50 M reads/minute).
=== Summary ===
Total reads processed: 26,038,229
Reads with adapters: 714,205 (2.7%) #读到的adapter数目
Reads written (passing filters): 26,038,229 (100.0%)
Total basepairs processed: 2,603,822,900 bp
Quality-trimmed: 82,577,636 bp (3.2%) #过滤了3.2%数据
Total written (filtered): 2,513,138,030 bp (96.5%) #保留了96.5%的数据
————————————————
版权声明：本文为CSDN博主「super_qun」的原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接及本声明。
原文链接：https://blog.csdn.net/weixin_44452187/java/article/details/87688276
- trim-galore并行处理时的生成文件顺序,参考文件：
trim-galore并行处理时的几个问题
五、 对水稻基因组建立索引（这步比较慢）
1. 构建索引 基本命令就是bowtie2-build 基因组序列.fa 索引名字
bowtie2-build all.fa rice
#all.fa: 水稻基因组序列（绝对路径）
#rice: 建立索引的前缀名
2. 索引后得到 6个文件 ，为.bt2索引文件的前缀
rice.1.bt2
rice.2.bt2
rice.3.bt2
rice.4.bt2
rice.rev.1.bt2
rice.rev.2.bt2
六、利用bowtie2比对
1. 参考：
Bowtie2使用方法与参数详细介绍
bowtie序列比对软件的使用
tophat2、hisat2、trimmomatic、cufflinks、samtools和stringtie
samtools的用法简介
SAMtools的使用说明
samtools用法详解
处理SAM、BAM你需要Samtools
2. 【方法】
nohup bowtie2 -p 8 \ #设置线程数
-x $index \ #索引路径
-1 ${fq1} -2 ${fq2} \ #双端测序fastq文件
-S SRR6795678.sam & #生成sam文件的文件名
3. bowtie相关参数
【必须参数】：
- -x <bt2-idx>
由bowtie2-build所生成的索引文件的前缀。首先 在当前目录搜寻，然后在环境变量BOWTIE2_INDEXES中制定的文件夹中搜寻。 - -1 <m1>
双末端测寻对应的文件1。可以为多个文件，并用逗号分开；多个文件必须和-2<m2>中制定的文件一一对应。比如:"-1flyA_1.fq,flyB_1.fq -2 flyA_2.fq,flyB_2.fq".测序文件中的reads的长度可以不一样。 - -2 <m2>
双末端测寻对应的文件2.-U <r>非双末端测寻对应的文件。可以为多个文件，并用逗号分开。测序文件中的reads的长度可以不一样。 - -S <hit>
所生成的SAM格式的文件前缀。默认是输入到标准输出。
4. 将sam 格式转为bam格式
【应用】
#写一个脚本
bowtie2 -p 8 -n 2 -x $index \
-1 ./SRR6795670_1_val_1.fq.gz \
-2 ./SRR6795670_2_val_2.fq.gz \
-S SRR6795670.sam #进行比对得到bam文件
samtools view -Sbh SRR6795669.sam > SRR6795669.bam
#得到的sam文件为bam文件
#-b 以BAM格式输出，可以用于samtools的后续分析
#-S 输入文件为SAM格式
#-u 以未压缩的BAM格式输出，可以节约时间，一般在管道执行时使用
#-h 在结果中包含头header
samtools sort SRR6795669.bam -@ 10 -o SRR6795669.sorted.bam #对bam文件进行排序
#-@ 线程数
samtools index -@ 10 SRR6795669.sorted.bam
5. 对bam文件进行QC
ls *sorted.bam |xargs -i samtools index {}
ls *sorted.bam | while read id ;do (samtools flagstat $id > $(basename $id ".sorted.bam").stat);done
#查看成功率：
cat *stat*|grep %
[thu@mn02 mapping]$ cat *stat*|grep %
51413901 + 0 mapped (71.24% : N/A)
45088678 + 0 properly paired (62.47% : N/A)
5349781 + 0 singletons (7.41% : N/A)
38939471 + 0 mapped (77.03% : N/A)
35879950 + 0 properly paired (70.98% : N/A)
2515027 + 0 singletons (4.98% : N/A)
67913650 + 0 mapped (93.74% : N/A)
67386764 + 0 properly paired (93.01% : N/A)
186796 + 0 singletons (0.26% : N/A)
67389488 + 0 mapped (94.44% : N/A)
66910562 + 0 properly paired (93.77% : N/A)
204310 + 0 singletons (0.29% : N/A)
七、samtoolsPCR重复的去除
#samtools rmdup 去重复
samtools rmdup input.sorted.bam rmdup.bam
#对rmdup.bam文件进行QC
ls *.rmdup.bam | while read id ;do (nohup samtools flagstat $id > $(basename $id ".bam").stat & );done
#查看QC质控结果：
cat *stat*|grep %
51413901 + 0 mapped (71.24% : N/A)
45088678 + 0 properly paired (62.47% : N/A)
5349781 + 0 singletons (7.41% : N/A)
38939471 + 0 mapped (77.03% : N/A)
35879950 + 0 properly paired (70.98% : N/A)
2515027 + 0 singletons (4.98% : N/A)
#去除重复前后比较（主要看文件大小）
ls -lh
八、利用macs2来callpeak
vim callpeaking #创建一个脚本文件名为callpeaking
1 #!bin/bash
2 #callpeak
3 #BSUB -J align145
4 #BSUB -n 10
5 #BSUB -R span[hosts=1]
6 #BSUB -o %J.out
7 #BSUB -e %J.err
8 #BSUB -q normal
9
10
11 control="/public/home/thu/chip_seq/sra/5.RemovePCRduplicated/input.rmdup.bam"
12 treatment="/public/home/thu/chip_seq/sra/5.RemovePCRduplicated/control_H3K9ac_Chipseq_rep2_rmdup.bam"
13 outname="control_H3K9ac_Chipseq_rep2_rmdup"
14
15 nohup macs2 callpeak -c ${control} -t ${treatment} -f BAM -g 400000000 -n ${outname} &
#-c 后面加上control组样本名
#-t 后面加上treatment组样本名
#-f {AUTO,BAM,SAM,BED,ELAND,ELANDMULTI,ELANDEXPORT,BOWTIE,BAMPE,BEDPE}
# 后面加上format 输出文件的格式
#-q QVALUE 后面加上人为设定的q值。默认为0.05
#-B, --bdg 加上这个选项则表示需要输出extended fragment pileup和local lambda tracks (two files)
#-n NAME 指定输出文件名
#-g GSIZE Effective genome size，可以是具体数字，也可以用简写：
# 'hs' for human (2.7e9), 'mm' for mouse (1.87e9), 'ce' for C. elegans (9e7)
# 'dm' for fruitfly (1.2e8), Default:hs
# 这里水稻的EGS为4.00e+08
# GES定义：A number of tools can accept an “effective genome size”. This is defined as the length of the “mappable” genome. There are two common alternative ways to calculate this:1. The number of non-N bases in the genome.2. The number of regions (of some size) in the genome that are uniquely mappable (possibly given some maximal edit distance).
#- 水稻基因组大小大约是4*10^8（总碱基数）。 一般要除去端粒和着丝粒等测序测不到的部位，但影响不大。
1. 【使用参数 】
-c: 对照组的输出结果
- control 或 mock(非特异性抗体，如IgG)组
control:input DNA，没有经过免疫共沉淀处理； - mock:
1）未使用抗体富集与蛋白结合的DNA片段
2）非特异性抗体，如IgG
-f: -t和-c提供文件的格式，可以是”ELAND”, “BED”, “ELANDMULTI”, “ELANDEXPORT”, “ELANDMULTIPET” (for pair-end tags), “SAM”, “BAM”, “BOWTIE”, “BAMPE” “BEDPE” 任意一个。如果不提供这项，就是自动检测选择。MACS2读入文件格式。
-g: 基因组大小， 默认提供了hs, mm, ce, dm选项， 不在其中的话，比如说拟南芥，就需要自己提供了。
- 有效基因组大小(可比对基因组大小)；基因组中有大量重复序列测序测不到，实际上可比对的基因组大小只有原基因组90% 或 70%；人类默认值是– 2.7e9（UCSC human hg18 assembly）
-n: 输出文件的前缀名
eg:NAME .会产生‘NAME_peaks.xls’, ‘NAME_negative_peaks.xls’, ‘NAME_peaks.bed’ , ‘NAME_summits.bed’, ‘NAME_model.r’四个文件
-B: 会保存更多的信息在bedGraph文件中，如fragment pileup, control lambda, -log10pvalue and -log10qvalue scores
-q: q值，也就是最小的PDR阈值， 默认是0.05。q值是根据p值利用BH计算，也就是多重试验矫正后的结果。
-p：这个是p值，指定p值后MACS2就不会用q值了。
-m: 和MFOLD有关，而MFOLD和MACS预构建模型有关，默认是5：50，MACS会先寻找100多个peak区构建模型，一般不用改，因为我们很大概率上不会懂。
2. macs2结果输出
control_H3K9ac_Chipseq_rep2_rmdup_model.r
control_H3K9ac_Chipseq_rep2_rmdup_peaks.narrowPeak
control_H3K9ac_Chipseq_rep2_rmdup_peaks.xls
control_H3K9ac_Chipseq_rep2_rmdup_summits.bed
3. 得到的 NAME_peaks.xls 文件
Chr2 17873714 17875118 1405 17874072 30 13.2471 5.04787 11.6898 control_H3K9ac_Chipseq_rep2_rmdup_peak_11539
Chr2 17880926 17882366 1441 17881114 47.02 30.3482 8.24866 28.4358 control_H3K9ac_Chipseq_rep2_rmdup_peak_11540
Chr2 17888892 17889834 943 17889240 55.13 37.9717 9.40114 35.9095 control_H3K9ac_Chipseq_rep2_rmdup_peak_11541
Chr2 17905088 17905946 859 17905509 69.73 54.316 12.1478 51.9143 control_H3K9ac_Chipseq_rep2_rmdup_peak_11542
虽然后缀名是xls，但实际上就是一个普通的文本文件。包含peak信息的tab分割的文件，前几行会显示callpeak时的命令。输出信息包含：
-
染色体号
-
peak起始位点
-
peak结束位点
-
peak区域长度
-
peak的峰值位点（summit position）
-
peak 峰值的属性（包括pileup峰高和可信度）（pileup height at peak summit, -log10(pvalue) for the peak summit）都是值越高越好
-
peak的富集倍数（相对于random Poisson distribution with local lambda）
4. 得到的 NAME_peaks.narrowPeak 文件
Chr1 28122534 28123365 control_H3K9ac_Chipseq_rep2_rmdup_peak_2602 555 . 12.3587 57.9451 55.5223 567
Chr1 28126628 28127930 control_H3K9ac_Chipseq_rep2_rmdup_peak_2603 233 . 7.42181 25.0983 23.3049 220
Chr1 28133217 28134651 control_H3K9ac_Chipseq_rep2_rmdup_peak_2604 211 . 6.82834 22.8667 21.1265 280
Chr1 28136385 28136908 control_H3K9ac_Chipseq_rep2_rmdup_peak_2605 93 . 4.4158 10.8192 9.32836 336
narrowPeak文件是BED6+4格式，可以上传到UCSC浏览。输出文件每列信息分别包含：
-
染色体号
-
peak起始位点
-
结束位点
-
peak name
-
int(-10*log10qvalue)可信度，值越大越好
-
正负链
-
fold change
-
-log10pvalue可信度，值越大越好
-
-log10qvalue可信度，值越大越好
-
relative summit position to peak start相对于起始位点来说peaks峰值的位置
bed格式的文件和其他Peak文件唯一不同的是——Peak文件中，由于以1为基，而bed文件是以0为基，所以peak的起始位置需要减1才是bed格式的文件。
5. 得到的 NAME_summits.bed 文件
Chr1 34159624 34159625 control_H3K9ac_Chipseq_rep2_rmdup_peak_3468 39.8115
Chr1 34162177 34162178 control_H3K9ac_Chipseq_rep2_rmdup_peak_3469 6.22312
Chr1 34165994 34165995 control_H3K9ac_Chipseq_rep2_rmdup_peak_3470 45.9992
Chr1 34179707 34179708 control_H3K9ac_Chipseq_rep2_rmdup_peak_3471 23.5895
BED格式的文件，包含peak的summits位置，第5列是-log10pvalue。如果想找motif，推荐使用此文件。
能够直接载入UCSC browser，用其他软件分析时需要去掉第一行。？
6. 得到的 NAME_model.r 文件
NAME_model.r能通过 $ Rscript NAME_model.r作图，得到是基于你提供数据的peak模型。
可以看到非常贴心的帮我们写好了脚本，直接运行即可出pdf图。
十、使用deeptools是进行可视化
1. 数据的连续密度图(在基因组浏览器上查看)
[1] bam文件转为bw文件
- 为什么要用bigWig(bw格式)? 主要是因为BAM文件比较大，直接用于展示时对服务器要求较大。
cd /public/home/thu/chip-seq/sra/5.RemovePCRduplicated
ls *.bam |xargs -i samtools index {} #进行转换前需要先建立index
cat >bam2bw.command #写一个脚本
#bin/bash
#callpeak
#BSUB -J align145
#BSUB -n 10
#BSUB -R span[hosts=1]
#BSUB -o %J.out
#BSUB -e %J.err
#BSUB -q normal
ls *.sorted.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.bw
done
ls *.rmdup.bam |while read id;do
bamCoverage -p 10 --normalizeUsing CPM -b $id -o ${id%%.*}.rmdup.bw
done
:wq
nohup bash bam2bw.command &
mkdir bwfile
mv *bw bwfile
cd bwfile
sz *.bw #将bw文件发送到本地
【2】载入IGV中进行查看：
- 用于可视化的基因组浏览器：IGV、JBrowse等
- IGV 中需要导入基因组参考文件和上面得到的bw文件：
genomes>load from file 选择对应的基因组参考文件
file>load from file 选择对应的bam文件 - 这里下载水稻基因序列参考文件all.con，网站 http://rice.plantbiology.msu.edu/，如下下载：
Rice Genome Annotation Project - 搜狗_51Record.gif
得到密度图：
image.png
2. 折线图/热图：所有转录本上peak区域reads的分布情况
【1】依赖：
- deeptools bamCoverage得到的bw文件
- 基因组注释文件bed文件
【2】原理：
- 为了统计全基因组范围的peak在基因特征的分布情况，用deeptools软件：也就是需要用到computeMatrix函数计算，用plotHeatmap函数以热图的方式对覆盖进行可视化，用plotProfile函数以折线图的方式展示覆盖情况。
【3】应用：
mkdir 7.rice_tss.bed && cd 7.rice_tss.bed
ln -s 5.RemovePCRduplicated/bwfile ./
#下载水稻基因组注释文件all.locus_brief_info.7.0（非bed格式）
nohup wget http://rice.plantbiology.msu.edu/pub/data/Eukaryotic_Projects/o_sativa/annotation_dbs/pseudomolecules/version_7.0/all.dir/all.locus_brief_info.7.0 &
#将下载的注释文件按顺序提取染色体号chrom、染色体起始位置chromStart、染色体终止位置chromEnd、基因名字四列，生成bed文件
cat all.locus_brief_info.7.0 | cut -f 1,4,5 > rice.1
cat all.locus_brief_info.7.0 | cut -f 2 > rice.2
paste rice1 rice2 > rice.bed
sed -i '1d' rice.bed
#创建computermatric脚本
vim computermatric
1 #!bin/bash
2 #computeMatric
3 #BSUB -J align145
4 #BSUB -n 10
5 #BSUB -R span[hosts=1]
6 #BSUB -o %J.out
7 #BSUB -e %J.err
8 #BSUB -q normal
9
10 bed=/public/home/thu/2.chip_seq/sra/7.rice_tss.bed/rice.bed
11 for id in bwfile/*bw;
12 do
13 file=$(basename $id )
14 sample=${file%%.bw}
15 echo $sample
16
17 computeMatrix reference-point --referencePoint TSS -p 15 \
18 -b 2000 -a 2000 \
19 -R $bed \
20 -S $id \
21 --skipZeros -o matrix1_${sample}_TSS_2K.gz \
22 --outFileSortedRegions regions1_${sample}_TSS_2K.bed
24 #plotHeatmap和plotProfile函数
25 plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.png
26 plotHeatmap -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Heatmap_2K.pdf --plotFileFormat pdf --dpi 720
27 plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.png
28 plotProfile -m matrix1_${sample}_TSS_2K.gz -out ${sample}_Profile_2K.pdf --plotFileFormat pdf --perGroup --dpi 720
29 done
#computeMatrix reference-point函数的参数
#--referencePoint {TSS,TES,center} The reference point for the plotting could be either the region start (TSS), the region end (TES) or the center of the region
#-b/--beforeRegionStartLength INT bp：区域前多少bp
#-a/--afterRegionStartLength INT bp：区域后多少bp
#-R/--regionsFileName File：指定bed文件
#-S/--scoreFileName File：指定bw文件
#--skipZeros：默认为False
#-o/-out/--outFileName OUTFILENAME：指定输出文件名
#--outFileSortedRegions BED file
#plotHeatmap和plotProfile函数参数
#-m/--matrixFile MATRIXFILE：指定matrix名字（computeMatrix的输出文件名）
#-o/-out/--outFileName OUTFILENAME：指定输出文件名，根据后缀自动决定图像的格式
#--plotFileFormat {"png", "eps", "pdf", "plotly"，"svg"}：指定输出文件格式
#--dpi DPI：指定dpi
#--perGroup：默认是画一个样本的所有区域群的图，使用这个选项后则按区域群画所有样本的情况
【4】结果如图：
image.png
【5】参考：
- peak注释信息揭秘 - 简书 https://www.jianshu.com/p/71028fd09d97
- 基因组注释文件格式 --（一）BED文件格式 https://mp.weixin.qq.com/s?src=11×tamp=1589444269&ver=2337&signature=lIzH-tC307nR24nyDIVhpyAt2kZ1RHm53erJdCyplJ5VVaMUcA7SVB-1cyt72nGL0hATEZjkHQC2ZrCtzVsyXuZYGpProA*8OzDZS6hCirKeWFZtPdvIR8we7X71oFkU&new=1
- deepTools 的使用（一） - 简书 https://www.jianshu.com/p/e7e2c65183fd
- ChIP-Seq数据挖掘系列-3: Motif 分析(3) - 利用ChIP-Seq结果在基因组区域中寻找富集的Motifs - 简书 https://www.jianshu.com/p/1602c2621c2b
网友评论