越哥cuttag分析代码留存

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 
最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
【社区内容提示】社区部分内容疑似由AI辅助生成,浏览时请结合常识与多方信息审慎甄别。
平台声明:文章内容(如有图片或视频亦包括在内)由作者上传并发布,文章内容仅代表作者本人观点,简书系信息发布平台,仅提供信息存储服务。

相关阅读更多精彩内容

友情链接更多精彩内容