1、trim_galore去接头及低质量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挂后台后可以用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文件夹下，会有val_1和val_2生成，同时有一个report
1694613373020.png
2、使用bowtie2比对
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、转格式
#sam转bam
samtools view -bS -F 0x04 H3K27ac-1.sam > H3K27ac-1.mapped.bam
#-F 0x04这个参数的意思就是转bam的时候去掉没比对上的那部分
#bam排序
samtools sort H3K27ac-1.mapped.bam -o H3K27ac-1.sort.bam
4、转bedGraph转bw
bedtools genomecov -bga -ibam H3K27ac-1.sort.bam > H3K27ac-1.bedGraph
#bedGraph排序后转bw
sort -k1,1 -k2,2n H3K27ac-1.bedGraph -o H3K27ac-1.sorted.bedGraph
#这里需要下载bedGraphToBigWig（wget下载，使用加绝对路径就行）
#还有基因组的hg38.fa需要建立index
samtools faidx hg38.fa
#会生成一个fai文件
bedGraphToBigWig H3K27ac-1.sorted.bedGraph /work/home/algroup01/genome/hg38_ucsc/hg38.fa.fai H3K27ac-1.bw
5、callpeak
callpeak这个需要根据不同的修饰不同的专利因子来调整参数
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指所有的转录因子
#对于H327ac是narrowpeak 默认窄峰
#对于PR这个转录因子，也是窄峰
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
网友评论