1、trim_galore去接頭及低質(zhì)量reads
nohup
trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /home/data/t170409/data/wwjtest/clean /home/data/t170409/data/oriadata/hg38_index/H3K27ac-2.R1.fq.gz /home/data/t170409/data/oriadata/hg38_index/H3K27ac-2.R2.fq.gz
&
使用nohup掛后臺(tái)后可以用jobs -l查看
jobs -l
[1] 1699330 Running nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /home/data/t170409/data/wwjtest/clean /home/data/t170409/data/oriadata/hg38_index/PR-6.R1.fq.gz /home/data/t170409/data/oriadata/hg38_index/PR-6.R2.fq.gz &
[2]- 1708775 Running nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /home/data/t170409/data/wwjtest/clean /home/data/t170409/data/oriadata/hg38_index/PR-4.R1.fq.gz /home/data/t170409/data/oriadata/hg38_index/PR-4.R2.fq.gz &
[3]+ 1710667 Running nohup trim_galore -q 25 --phred33 --length 35 -e 0.1 --stringency 4 --paired -o /home/data/t170409/data/wwjtest/clean /home/data/t170409/data/oriadata/hg38_index/PR-3.R1.fq.gz /home/data/t170409/data/oriadata/hg38_index/PR-3.R2.fq.gz &
跑完以后在clean文件夾下,會(huì)有val_1和val_2生成听盖,同時(shí)有一個(gè)report
2杨拐、使用bowtie2比對(duì)
nohup bowtie2 -p 8 --very-sensitive-local --no-unal --no-mixed --no-discordant --phred33 -I 10 -X 700 -x /home/data/t170409/data/oriadata/hg38_index/hg38 -1 /home/data/t170409/data/wwjtest/clean/H3K27ac-2.R1_val_1.fq.gz -2 /home/data/t170409/data/wwjtest/clean/H3K27ac-2.R2_val_2.fq.gz -S /home/data/t170409/data/wwjtest/align/H3K27ac-2.sam > /home/data/t170409/data/wwjtest/align/H3K27ac-2.log &
3棕所、轉(zhuǎn)格式
#sam轉(zhuǎn)bam
samtools view -bS -F 0x04 H3K27ac-1.sam > H3K27ac-1.mapped.bam
#-F 0x04這個(gè)參數(shù)的意思就是轉(zhuǎn)bam的時(shí)候去掉沒(méi)比對(duì)上的那部分
#bam排序
samtools sort H3K27ac-1.mapped.bam -o H3K27ac-1.sort.bam
4衫生、轉(zhuǎn)bedGraph轉(zhuǎn)bw
bedtools genomecov -bga -ibam H3K27ac-1.sort.bam > H3K27ac-1.bedGraph
#bedGraph排序后轉(zhuǎn)bw
sort -k1,1 -k2,2n H3K27ac-1.bedGraph -o H3K27ac-1.sorted.bedGraph
#這里需要下載bedGraphToBigWig(wget下載裳瘪,使用加絕對(duì)路徑就行)
#還有基因組的hg38.fa需要建立index
samtools faidx hg38.fa
#會(huì)生成一個(gè)fai文件
bedGraphToBigWig H3K27ac-1.sorted.bedGraph /work/home/algroup01/genome/hg38_ucsc/hg38.fa.fai H3K27ac-1.bw
5、callpeak
callpeak這個(gè)需要根據(jù)不同的修飾不同的專(zhuān)利因子來(lái)調(diào)整參數(shù)
H3F3A','H3K4me1','H3K9me1','H4K20me1','H3K27me3','H3K79me2','H3K9me2','H3K36me','H3K79me3','H3K9me3'這些都是broad寬峰罪针,'TF','H2AFZ','H3K27ac','H3K4me3','H3ac','H3K9ac','H3K4me2','H2A.Z','RNA_pol_II','IgG'這些都是narrow窄峰,TF指所有的轉(zhuǎn)錄因子
#對(duì)于H327ac是narrowpeak 默認(rèn)窄峰
#對(duì)于PR這個(gè)轉(zhuǎn)錄因子黄伊,也是窄峰
macs2 callpeak -f BAM -t H3K27ac-2.sort.bam -n "H3K27ac-2" -g hs --nomodel --keep-dup all --outdir /work/home/algroup01/username/wangwenjing/wangyue_cuttag/H3K27ac/callpeak