這里是佳奧！我們繼續(xù)clean數據的比對。

1 使用bowtie2進行比對

用bowtie2進行比對和統(tǒng)計比對率, 需要提前下載參考基因組然后使用命令構建索引，或者直接就下載索引文件：

下載小鼠參考基因組的索引和注釋文件, 這里用常用的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功能，避免麻煩。
## http://www.itdecent.cn/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

了解相關參數

http://www.itdecent.cn/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ū)域的局部偏差/λ

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的一半。
• 示例：

1. 想找富集剪切位點，如DNAse-seq，所有5'端的序列reads應該從兩個方向延伸，如果想設置移動的窗口是200bp，參數設置如下：
   --nomodel --shift -100 --extsize 200
2. 對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計算重復情況以及可視化的內容。

我們下一篇再見！