这里是佳奥！我们继续clean数据的比对。
1 使用bowtie2进行比对
用bowtie2进行比对和统计比对率, 需要提前下载参考基因组然后使用命令构建索引，或者直接就下载索引文件：
下载小鼠参考基因组的索引和注释文件, 这里用常用的mm10
# 索引大小为3.2GB，不建议自己下载基因组构建，可以直接下载索引文件，代码如下：
mkdir referece && cd reference
wget -4 -q ftp://ftp.ccb.jhu.edu/pub/data/bowtie2_indexes/mm10.zip
unzip mm10.zip
如果要判断属于哪一个物种：选取一段序列blast
https://genome.ucsc.edu/cgi-bin/hgBlat?command=start
##对bam文件要进行以下处理
1 去除PCR重复
2 去除低质量reads
3 去除未比对到同一染色体
4 去除线粒体
小样本测试
zcat ../clean/2-cell-1_1_val_1.fq.gz |head -10000 > test1.fq
zcat ../clean/2-cell-1_2_val_2.fq.gz |head -10000 > test2.fq
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq -S test.sam
bowtie2 -x /public/reference/index/bowtie/mm10 -1 test1.fq -2 test2.fq |samtools sort -@ 5 -O bam -o test.bam -
## 建议抛弃 samtools markdup功能，避免麻烦。
## https://www.jianshu.com/p/1e6189f641db
samtools markdup -r test.bam test.samtools.rmdup.bam
## 把报错信息在谷歌搜索后，在两个网页上找到了答案。
https://github.com/samtools/samtools/issues/765
https://www.biostars.org/p/288496/
## gatk 可以在GitHub下载
/public/biosoft/GATK/gatk-4.0.3.0/gatk MarkDuplicates \
-I test.bam -O test.picard.rmdup.bam --REMOVE_SEQUENCING_DUPLICATES true -M test.log
### picards 被包装在GATK里面：
### sambamba 文档： http://lomereiter.github.io/sambamba/docs/sambamba-markdup.html
conda install -y sambamba
sambamba --help
sambamba markdup --help
sambamba markdup -r test.bam test.sambamba.rmdup.bam
samtools flagstat test.sambamba.rmdup.bam
samtools flagstat test.bam
## 接下来只保留两条reads要比对到同一条染色体(Proper paired) ，还有高质量的比对结果(Mapping quality>=30)
## 顺便过滤 线粒体reads
samtools view -f 2 -q 30 test.sambamba.rmdup.bam |grep -v chrM|wc
samtools view -f 2 -q 30 test.sambamba.rmdup.bam |wc
samtools view -h -f 2 -q 30 test.sambamba.rmdup.bam |grep -v chrM| samtools sort -O bam -@ 5 -o - > test.last.bam
bedtools bamtobed -i test.last.bam > test.bed
ls *.bam |xargs -i samtools index {}
正式比对脚本
ls /home/kaoku/project/atac/clean/*_1.fq.gz > 1
ls /home/kaoku/project/atac/clean/*_2.fq.gz > 2
ls /home/kaoku/project/atac/clean/*_2.fq.gz | cut -d"/" -f 7 | cut -d"_" -f 1 > 0
paste 0 1 2 > config.clean ##供mapping使用的配置文件
##相对目录需要理解
cd ~/project/atac/align
##一定要搞清楚自己的bowtie2软件安装在哪里，以及自己的索引文件在什么地方
bowtie2_index=/home/kaoku/refer/mm10/mm10
cat config.clean |while read id;
do echo $id
arr=($id)
fq2=${arr[2]}
fq1=${arr[1]}
sample=${arr[0]}
##比对过程15分钟一个样本
bowtie2 -p 5 --very-sensitive -X 2000 -x $bowtie2_index -1 $fq1 -2 $fq2 |samtools sort -O bam -@ 5 -o - > ${sample}.raw.bam
samtools index ${sample}.raw.bam
bedtools bamtobed -i ${sample}.raw.bam > ${sample}.raw.bed
samtools flagstat ${sample}.raw.bam > ${sample}.raw.stat
# https://github.com/biod/sambamba/issues/177
sambamba markdup --overflow-list-size 600000 --tmpdir='./' -r ${sample}.raw.bam ${sample}.rmdup.bam
samtools index ${sample}.rmdup.bam
## ref:https://www.biostars.org/p/170294/
## Calculate %mtDNA:
## mtReads=$(samtools idxstats ${sample}.rmdup.bam | grep 'chrM' | cut -f 3)
## totalReads=$(samtools idxstats ${sample}.rmdup.bam | awk '{SUM += $3} END {print SUM}')
## echo '==> mtDNA Content:' $(bc <<< "scale=2;100*$mtReads/$totalReads")'%'
samtools flagstat ${sample}.rmdup.bam > ${sample}.rmdup.stat
samtools view -h -f 2 -q 30 ${sample}.rmdup.bam |grep -v chrM |samtools sort -O bam -@ 5 -o - > ${sample}.last.bam
samtools index ${sample}.last.bam
samtools flagstat ${sample}.last.bam > ${sample}.last.stat
bedtools bamtobed -i ${sample}.last.bam > ${sample}.bed
done
其中bowtie2比对加入了-X 2000 参数，是最大插入片段，宽泛的插入片段范围(10-1000bp)
上述脚本的步骤都可以拆分运行，比如bam文件构建index或者转为bed的：
ls *.last.bam|xargs -i samtools index {}
ls *.last.bam|while read id;do (bedtools bamtobed -i $id >${id%%.*}.bed) ;done
ls *.raw.bam|while read id;do (nohup bedtools bamtobed -i $id >${id%%.*}.raw.bed & ) ;done
##宽泛的插入片段范围(10-1000bp)
两条reads要比对到同一条染色体(Proper paired)
高质量的比对结果(Mapping quality> =30)
常用软件: bowtie2、bwa
命令实例:
- Bowtie2 -x index -X 1000 -1 read1.fq -2 read2.fq -S output.sam
- Samtools view -f 2 -q 30 -o output.filter.bam output.sam
##PCR重复，前后质控
线粒体、叶绿体等细胞器来源的数据
需要软件:picard, samtools
java -jar picard.jar MarkDuplicates I=output.filter.bam O=output.dedup.bam
M=duplication.log
查看线粒体数据的比例,去除线粒体污染
- samtools view -h output.dedup.bam | grep -v 'chrM' | samtools view -bS -o output.final.bam
过滤后，比对得到文件
-rw-r--r-- 1 root root 237M 12月 30 12:34 2-ceLL-1.bed
-rw-r--r-- 1 root root 488M 12月 30 12:33 2-ceLL-1.last.bam
-rw-r--r-- 1 root root 2.9M 12月 30 12:33 2-ceLL-1.last.bam.bai
-rw-r--r-- 1 root root 3.7G 12月 30 12:19 2-ceLL-1.raw.bam
-rw-r--r-- 1 root root 2.9M 12月 30 12:20 2-ceLL-1.raw.bam.bai
-rw-r--r-- 1 root root 2.5G 12月 30 12:21 2-ceLL-1.raw.bed
-rw-r--r-- 1 root root 778M 12月 30 12:33 2-ceLL-1.rmdup.bam
-rw-r--r-- 1 root root 2.8M 12月 30 12:33 2-ceLL-1.rmdup.bam.bai
##查看.bam文件
samtools idxstats 2-ce11-2.raw.bam | grep -v "_"
chr1 195471971 15013607 13329
chr2 182113224 2104122 7406
chr3 160039680 774572 2919
chr4 156508116 771827 3531
chr5 151834684 839608 2962
chr6 149736546 833724 3460
chr7 145441459 741301 3202
chr8 129401213 645650 2711
chr9 124595110 741435 2838
chr10 130694993 731066 2238
chr11 122082543 788663 3755
chr12 120129022 981501 5700
chr13 120421639 663876 2267
chr14 124902244 662372 2918
chr15 104043685 518410 2352
chr16 98207768 524718 1879
chr17 94987271 517792 2470
chr18 90702639 476631 2106
chr19 61431566 372339 1658
chrX 171031299 737566 1774
chrY 91744698 58961 583
chrM 16299 65542629 214917 ##线粒体
##查看两次过滤去除情况
samtools idxstats 2-ce11-2.raw.bam | grep -v "_" > 1
samtools idxstats 2-ce11-2.rmdup.bam | grep -v "_" > 2
samtools idxstats 2-ce11-2.last.bam | grep -v "_" > 3
paste 1 2 3 | cut -f 1,3,7,11
chr1 15013607 1326747 565438
chr2 2104122 692987 554436
chr3 774572 465602 406922
chr4 771827 463522 395260
chr5 839608 496581 424198
chr6 833724 451431 395648
chr7 741301 459578 380588
chr8 645650 398156 339770
chr9 741435 446255 377214
chr10 731066 436419 389050
chr11 788663 480939 434394
chr12 981501 403747 326044
chr13 663876 398679 341710
chr14 662372 384470 314884
chr15 518410 314024 281836
chr16 524718 316711 283482
chr17 517792 317753 268980
chr18 476631 286928 257332
chr19 372339 225302 204708
chrX 737566 454629 336846
chrY 58961 48319 3136
chrM 65542629 2854549 0 ##线粒体
* 0 0
2 使用macs2找peaks
了解相关参数
https://www.jianshu.com/p/21e8c51fca23
输入文件参数：
-
-t:实验组，IP的数据文件 -
c: 对照组 -
f：指定输入文件的格式，默认是自动检测输入数据是什么格式，支持bam,sam,bed等 - g：有效基因组大小，由于基因组序列的重复性，基因组实际可以mapping的大小小于原始的基因组。这个参数要根据实际物种计算基因组的有效大小。软件里也给出了几个默认的-g 值：hs -- 2.7e9表示人类的基因组有效大小(UCSC human hg18 assembly).
- hs: 2.7e9
- mm: 1.87e9
- ce: 9e7
- dm: 1.2e8
输出文件参数：
-
--outdir:输出结果的存储路径
--n:输出文件名的前缀 -
-B/--bdg:输出bedgraph格式的文件，输出文件以NAME+'_treat_pileup.bdg' for treatment data, NAME+'_control_lambda.bdg' for local lambda values from control显示。
peak calling 参数
-
-q/--qvalue和-p/--pvalue
q value默认值是0.05，与pvalue不能同时使用。 -
--broad
peak有narrow peak和broad peak, 设置时可以call broad peak 的结果文件。 -
--broad-cutoff
和pvalue、以及qvalue相似 -
--nolambda: 不要考虑在峰值候选区域的局部偏差/λ
q值与峰宽有一定的联系。理想情况下，如果放宽阈值，您将简单地获得更多的峰值，但是使用MACS2放松阈值也会导致更宽的峰值。
Shift 模型参数：
-
--nomodel
这个参数和extsize、shift是配套使用的，有这个参数才可以设置extsize和shift。 -
--extsize
当设置了nomodel时，MACS会用--extsize这个参数从5'->3'方向扩展reads修复fragments。比如说你的转录因子结合范围200bp，就设置这个参数是200。 -
--shift
当设置了--nomodel，MACS用这个参数从5' 端移动剪切，然后用--extsize延伸，如果--shift是负值表示从3'端方向移动。建议ChIP-seq数据集这个值保持默认值为0，对于检测富集剪切位点如DNAsel数据集设置为EXTSIZE的一半。 - 示例：
- 想找富集剪切位点，如DNAse-seq，所有5'端的序列reads应该从两个方向延伸，如果想设置移动的窗口是200bp，参数设置如下：
--nomodel --shift -100 --extsize 200 - 对nucleosome-seq数据，用核小体大小的一半进行extsize,所以参数设置如下：
--nomodel --shift 37 --extsize 73
-
--call-summits
MACS利用此参数重新分析信号谱，解析每个peak中包含的subpeak。对相似的结合图谱，推荐使用此参数，当使用此参数时，输出的subpeak会有相同的peak边界，不同的绩点和peak summit poisitions.
批量处理脚本
##对一个样本处理
macs2 callpeak -t 2-cell-1.bed -g mm --nomodel --shift -100 --extsize 200 -n 2-cell-1 --outdir ../peaks/
ls *.bed | while read id ;
do (macs2 callpeak -t $id -g mm --nomodel --shift -100 --extsize 200 -n ${id%%.*} --outdir ../peaks/) ;
done
得到的文件
2-ce11-2_peaks.narrowPeak 2-ce11-4_peaks.narrowPeak 2-ce11-5_peaks.narrowPeak 2-ceLL-1_peaks.narrowPeak
2-ce11-2_peaks.xls 2-ce11-4_peaks.xls 2-ce11-5_peaks.xls 2-ceLL-1_peaks.xls
2-ce11-2_summits.bed 2-ce11-4_summits.bed 2-ce11-5_summits.bed 2-ceLL-1_summits.bed
根据.narrowPeak结果进到IGV查看
$ cat 2-ce11-2_peaks.narrowPeak | head
chr1 3175662 3176103 2-ce11-2_peak_1 86 . 7.76699 13.34822 8.64135 170
chr1 3445564 3445764 2-ce11-2_peak_2 44 . 5.69363 8.07027 4.47755 93
chr1 3930951 3931169 2-ce11-2_peak_3 35 . 4.78261 6.86304 3.57156 82
chr1 3953679 3953879 2-ce11-2_peak_4 27 . 4.36893 5.74423 2.72116 61
chr1 3958963 3959552 2-ce11-2_peak_5 128 . 8.07453 18.52744 12.84147 195
chr1 4448079 4448353 2-ce11-2_peak_6 63 . 6.81818 10.41459 6.31841 129
chr1 4699234 4699987 2-ce11-2_peak_7 59 . 5.71429 9.96971 5.96806 371
chr1 4712462 4712699 2-ce11-2_peak_8 32 . 4.71698 6.50781 3.28730 56
chr1 4718657 4718895 2-ce11-2_peak_9 41 . 5.12821 7.69019 4.19980 118
chr1 4760132 4760526 2-ce11-2_peak_10 35 . 4.74138 6.79799 3.51876 60
QQ截图20221231152406.png
可以看到有峰。
下一步便是计算插入片段长度，FRiP值，IDR计算重复情况以及可视化的内容。
我们下一篇再见！
