一、Chip-seq分析流程
二、数据、参考基因组、所需软件下载
1、参考基因组下载:
以二穗短柄草为例:
https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Bdistachyon
https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Bdistachyon#
2、所需软件下载:
①质控所需软件:fastqc
②过滤所需软件:trim_galore、trimmomatic
③去PCR重复所需软件:fastuniq
④比对所需软件:bowtie2、samtools、bedtools
⑤call peak所需软件:SICER
⑥绘图所需软件:R或者Rstudio
⑦可视化工具:IGV
前五个所需软件均可以使用conda安装,见上一篇文章《Conda 安装软件万能链接》:Conda安装软件万能链接
三、各步运行脚本
1、fastqc
fastqc x.fastq.gz 注:x.fastq.gz是下机数据,也就是raw data
结果是网页文件,可以下载到本地用浏览器打开看结果,后面会更新详细的结果说明文章
2、trim_galore
trim_galore --paired --phred33 --gzip x_1.fq.gz x_2.fq.gz
注:现在一般都是双端测序,所以有一个样品会得到两个测序文件
3、trimmomatic
java -jar trimmomatic-0.38.jar PE \
-threads 10 \
-phred33 \
X_R1.fastq.gz X_R2.fastq.gz \
X_paired_R1.fq.gz X_unpaired_R1.fq.gz X_paired_R2.fq.gz
X_unpaired_
R2.fq.gz \
ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:8:true \
SLIDINGWINDOW:5:20 \
LEADING:3 \
TRAILING:3 \
HEADCROP:13 \
MINLEN:36
4、fastuniq
fastuniq -i xy.txt -o x_uniq_R1.fq -p x_uniq_R2.fq
xy.txt文件中包含输入序列的名字,一行一个名字,格式如:
x_R1_1.fastq
y_R1_2.fastq
5、bowtie2
①建索引:bowtie2-build -f z.fa z_index
②比对:bowtie2 -p 10 -x z_index -1 x_paired_R1.fq -2 x_paired_R2.fq -S x.sam
6、samtools
转化格式:samtools view -bS -q 20 x.sam >x.bam
查看:samtools flagstat x.bam
排序:samtools sort x.bam x_sort.bam
建索引:samtools index x_sort.bam
7、bedtools
bedtools bamtobed -i x.bam >x.bed
8、SICER
sh SICER.sh x_exp_sort.bed x_input_sort.bed sicer_results species 1 200 150 0.8 200 0.01
9、划窗口
samtools faidx species.fa
awk '{print $1"\t"$2}' species.fa.fai >species.txt
bedtools makewindows -g species.txt -w chk >species_chk.bed
注:-w chk是指窗口大小,可以根据需要自己制定,如1kb窗口的话则-w chk替换成-w 1000
10、获得绘图数据
bedtools coverage -a species_chk.bed -b x_sort.bam >x_maping.txt
四、R包绘图
install.packages("ggplot2")
library(ggplot2)
read.table("D:x_maping.txt")
data<-read.table("D:x_maping.txt",header=T,sep="\t",na.strings="")
df <- ggplot(data)
df+aes(x=begin, y = mapping )+geom_bar(fill="red", stat = "identity",position = "dodge")+facet_grid(chr~.,space="free_x", scales="free_x")+ylim(0,50))
还有很多不会的地方,欢迎大家批评指正!