一法焰、Chip-seq分析流程
二址愿、數(shù)據(jù)混移、參考基因組祠墅、所需軟件下載
1、參考基因組下載:
以二穗短柄草為例:
https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Bdistachyon
https://genome.jgi.doe.gov/portal/pages/dynamicOrganismDownload.jsf?organism=Bdistachyon#
2沫屡、所需軟件下載:
①質(zhì)控所需軟件:fastqc
②過濾所需軟件:trim_galore、trimmomatic
③去PCR重復(fù)所需軟件:fastuniq
④比對所需軟件:bowtie2撮珠、samtools沮脖、bedtools
⑤call peak所需軟件:SICER
⑥繪圖所需軟件:R或者Rstudio
⑦可視化工具:IGV
前五個所需軟件均可以使用conda安裝,見上一篇文章《Conda 安裝軟件萬能鏈接》:Conda安裝軟件萬能鏈接
三芯急、各步運(yùn)行腳本
1勺届、fastqc
fastqc x.fastq.gz? 注:x.fastq.gz是下機(jī)數(shù)據(jù),也就是raw data
結(jié)果是網(wǎng)頁文件娶耍,可以下載到本地用瀏覽器打開看結(jié)果免姿,后面會更新詳細(xì)的結(jié)果說明文章
2、trim_galore
trim_galore --paired --phred33 --gzip x_1.fq.gz x_2.fq.gz
注:現(xiàn)在一般都是雙端測序榕酒,所以有一個樣品會得到兩個測序文件
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
轉(zhuǎn)化格式: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是指窗口大小肢础,可以根據(jù)需要自己制定还栓,如1kb窗口的話則-w chk替換成-w 1000
10、獲得繪圖數(shù)據(jù)
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))
還有很多不會的地方剩盒,歡迎大家批評指正!