DNA-seq

超详Chip-seq分析流程

2021-03-20  本文已影响0人  想要学好生信的小白

一、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安装软件万能链接

R安装链接:https://www.r-project.org/


Rstudio下载地址:https://www.rstudio.com/products/rstudio/download/


三、各步运行脚本

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))

还有很多不会的地方,欢迎大家批评指正!

上一篇下一篇

猜你喜欢

热点阅读