聲明：因為本人使用過程中遇到了很多原來人碰得到的問題，放防止以后遺忘，所以做一下記錄，因為前總結的很到位，所以直接引用到自己日常記錄，供日后查閱。
主要引用為：【原創(chuàng)文章】辣椒BSA分析對數據進行質控、去重、索引等預處理的幾個姿勢 （https://www.baishujun.com/archives/7191.html）！
去除PCR擴增重復序列
由于構建文庫的時候進行了pcr擴增，然后進行測序，這么測序出來的數據存在很多重復序列，在這里需要去除，然后建立索引，那么這個預處理工作就結束了。
遇到的問題最多的也是在這里。
按照之前的資料，對于PE測序數據一般用picard，效果好，SE可以用samtools或picard去重。
所以我首先用picard的MarkDuplicates進行了去重，遇到N個錯誤！
1，臨時目錄不夠
錯誤提示：
java運行問題，原因/tmp空間太小造成的，有兩個解決辦法。
（1），給picard設置特定的tmp目錄 (原作者推薦）
picard MarkDuplicates I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt TMP_DIR=/www/tmp
本人項目中也是用特定目錄解決， 在picard運行命令之前，加一行命令：
export _JAVA_OPTIONS=-Djava.io.tmpdir=./tmp
picard -Xmx32g MarkDuplicates I=S99.new1.bam O=S99.rmdup.bam CREATE_INDEX=ture REMOVE_DUPLICATES=ture M=S99.marked_dup_metrics.txt
（2），為了一勞永逸的解決這個問題，我直接掛載了一個1t的分區(qū)給/tmp（原作者）
#####此處還有一點項目中的查詢結果針對 java.lang.OutOfMemoryError: GC overhead limit exceeded
如果在垃圾收集中花費太多的時間太少，GC會拋出這個異常。GC上的CPU時間占用了98％，只有不到2％的堆被恢復。
此功能旨在防止應用程序長時間運行，而由于堆太小，進行很少或沒有進度。
您可以使用命令行選項關閉此功能? -XX:-UseGCOverheadLimit （我在使用beagle嘗試使用過這個但是沒有解決問題，其他原因）。
2，java堆大小設置過低
錯誤提示：
面對“Java heap space”的問題，主要是由于我的測序數據大，占用內存大，但是java程序默認的堆大小分配不夠，所以出現這個問題，我們可以給當前運行的java程序設定一個內存值，從而覆蓋默認設置。
picard -Xmx40g MarkDuplicates I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt （Xmx40g，設置JVM最大可用內存為40G.）
3，“Insert size out of range”錯誤
關于這個錯誤，從字面意思上來看是插入大小超過了范圍，并且我發(fā)現基本都是在檢測重復完成之后才出現，這個錯誤應當和程序讀取reads的大小或者java寫入磁盤文件大小相關方面的限制有關，但是沒有找到具體原因，從錯誤提示也可以看出，這是在確認sam數據時出現的問題，那么我們可以采用寬松一點的審計策略，就可以繞過這個“Insert size out of range”錯誤。
picard -Xmx40g MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt
我用的是conda環(huán)境，所以上面的命令是我運行picard進行序列去重成功的命令。
如果不想使用conda的picard，那么我們也可以直接調用java包，在picard命令中已經包含了當前picard的jar包的地址，我們可以修改命令如下：
java -Xmx40g -jar /www/soft/miniconda2/envs/pepper/share/picard-2.20.3-0/picard.jar MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt
使用了VALIDATION_STRINGENCY=LENIENT參數以后
我們可以看到大量的提示
?其它去重軟件
經過這一輪之后，整理發(fā)現，還有好幾個很優(yōu)秀的去重軟件，
1，Sambamba
地址：http://lomereiter.github.io/sambamba/
原文描述如下：
Marks (by default) or removes duplicate reads. For determining whether a read is a duplicate or not, the same criteria as in Picard are used.
這是我推薦的替代picard的軟件
她篩選reads是否重復的標準和picard是一樣的（For determining whether a read is a duplicate or not, the same criteria as in Picard are used）。
關鍵是這個軟件支持多線程，速度就快多了，她也支持是標記或者刪除重復序列，并且自動建立bai索引文件。
該軟件還有其它的sort，index，merge等功能，詳細地址為，http://lomereiter.github.io/sambamba/docs/sambamba-markdup.html#OPTIONS，
需要代理才能打開。
sambamba markdup -p -t 16 S347.merged.sorted.viewed.bam S347.sorted.viewed.markdupba.bam
-p顯示過程，-t線程數16.
對于一個46.70GB的bam文件進行去重，picard耗時78.11分鐘，而sambamba僅僅用時33分鐘，快了一半吧。
但是我用picard去重標記重復后的文件為48.23G，用Sambamba去重后的文件為46.27g，這是否預示著picard標記了更多重復？
運行sambamsa軟件，我遇到一個Too many open files問題
sambamba-markdup: Cannot open file `/tmp/sambamba-pid87983-markdup-xrmt/PairedEndsInfozkeh29' in mode `w+' (Too many open files)
這是由于我們的系統(tǒng)限制導致的，解決方法有兩個。
（1），臨時辦法
首先
運行命令
ulimit -n
查看線程限制，默認一般為1024，我們修改為10240.
運行命令
ulimit -n 10240
（2），上面的方法重啟后就失效了，如果想重啟后生效，需要修改/etc/security/limits.conf文件，root權限。
在文件末尾添加
* soft nofile 10240
* hard nofile 10240
*代表所有該系統(tǒng)的用戶
這個10240值不能大于/proc/sys/fs/nr_open文件中的數值，否則將無法正常登錄系統(tǒng)。
我的nr_open值為1048576
2，samtools
samtools也有標記或者刪除重復序列的功能，命令分別為
samtools rmdup
刪除重復序列，但是軟件還是推薦標記重復序列
samtools markdup -@ 8 S347.merged.sorted.viewed.bam S347.sorted.viewed.markdupba.bam
通過資料查詢，發(fā)現這個工具對于bam有一定要求：bam文件以染色體位置排序，MC tag需要用fixmate -m來提供，fixmate -m操作的bam文件必須是利用read name來排序的。
我放棄了，沒有測試。
3，samblaster
地址：https://github.com/GregoryFaust/samblaster
命令示例：
samblaster -i samp.disc.sam -o samp.split.sam
沒有多線程，操作的是sam文件，沒有測試。
4，biobambam2
地址：
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4075596/
https://github.com/gt1/biobambam2
bammarkduplicates功能可以去重，沒測試
5，bamUtil
https://github.com/statgen/bamUtil
有興趣的試試吧
生成索引
1，用Sambamba去重時，會自動生成bai索引文件。
但是當我使用命令
sambamba index -p S347.sorted.viewed.markdup.bam
卻顯示錯誤
Segmentation fault (core dumped)
2，picard MarkDuplicates去重時加入CREATE_INDEX=true參數，也可以建立索引。
單獨使用picard BuildBamIndex命令
picard BuildBamIndex -Xmx64g I=S352.sorted.viewed.markdup.bam
對于picard MarkDuplicates生成的bam提示“非bam”錯誤
picard.sam.BuildBamIndex done. Elapsed time: 0.00 minutes.
Runtime.totalMemory()=2147483648
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread “main” htsjdk.samtools.SAMException: Input file must be bam file, not sam file.
這個真是搞不懂了
用flagstat檢查了一下，好像bam文件正常啊
samtools flagstat -@ 30 S346.sorted.viewed.markdup.bam >S346.sorted.markdup.stat
對于sambamba markdup生成的bam，不說不是bam文件了，但是提示“ Insert size out of range”錯誤
picard.sam.BuildBamIndex done. Elapsed t ime: 29.70 minutes.
Runtime.totalMemory()=5200936960
To get help, see http://broadinstitute.github.io/picard/index.html#Gett ingHelp
Exception in thread “main” htsjdk.samtools.SAMFormatException: SAM vali dation error: ERROR: Record 704691907, Read name A00174:10:HCCJKDSXX:4: 2426:2564:33301, Insert size out of range
at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.ja va:455)
3，使用samtools index建立索引
samtools index S347.sorted.viewed.markdup.bam
卻提示錯誤
Region 536874585..536874735 cannot be stored in a bai index. Try using a csi index with min_shift = 14, n_lvls >= 6
samtools index: failed to create index for “S347.sorted.viewed.markdup.bam”: Numerical result out of range
查詢了一下資料，原始BAM索引格式（BAI）具有固定的bin大小系統(tǒng)，并且最大參考序列長度為512Mbp，使其不適合大染色體。
辣椒的參考基因組cm334，zunla超過3個G，遠遠大于536,870,912，bai格式索引不適合辣椒基因組分析！
辣椒只能使用csi格式索引文件，然后使用bcftools來call snp了
samtools index -c S347.sorted.viewed.markdup.bam
把變異檢測前對數據進行預處理的命令總結一下：
#數據指控
fastp -q 15 -u 40 -M 0 -n 5 -l 80 -w 16 -i? S347_jikangda-A_Green-beifen-1_AHCCJKDSXX_S347_L004_R1_001.fastq.gz -I? S347_jikangda-A_Green-beifen-1_AHCCJKDSXX_S347_L004_R2_001.fastq.gz -o S347.clean_R1.fastq.gz -O S347.clean_R2.fastq.gz -h report/S347-report.html
#建立參考基因組索引
bwa index Capsicum.annuum.L_Zunla-1_Release_2.0.fasta -p zunla
#與參考基因組比對并轉換成bam文件
bwa mem -k 32 -M -t 16 -R '@RG\tID:S347\tPL:illumina\tLB:S347\tSM:S347'? zunla? S347.clean_R1.fastq.gz? S347.clean_R2.fastq.gz | samtools view -S -b - > S347.bam
#對bam文件排序
samtools sort -@ 10 -m 80G -O bam -o S347.sorted.viewed.bam S347.bam
#merge兩個lanes
samtools merge S347.merged.sorted.viewed.bam S347.sorted.viewed.bam S411.sorted.viewed.bam （此處samtools 其實可以多線程）
#去除PCR重復序列
picard -Xmx40g? MarkDuplicates VALIDATION_STRINGENCY=LENIENT I=S347.merged.sorted.viewed.bam O=S347.sorted.viewed.markdup.bam M=S347.sorted.viewed.markdup_metrics.txt
#構建bam文件的csi索引
samtools index -c S347.sorted.viewed.markdup.bam
