BSA分析查找极端差异表现基因的方法，与参考基因组进行比对当前样本的差异变化，首先需要对测序回来的数据进行数据质控、与参考基因组进行比对、然后检测编译，最后绘图。
一，数据质控
我用的fastp，多线程，和html报告
fastp -q 15 -u 40 -M 0 -n 5 -l 80 -w 16 -i S413_gencodon-A_10_AHFTMFDSXX_S3_L003_R1_001.fastq.gz -I S413_gencodon-A_10_AHFTMFDSXX_S3_L003_R2_001.fastq.gz -o S413.clean_R1.fastq.gz -O S413.clean_R2.fastq.gz -h report/S413-report.html
备注：输入两个raw reads文件，输出两个clean reads文件，gzip压缩。程序运行线程数16，期望碱基质量值设定为q15（-q 15）并去除高于40%低质量碱基（-u 40）数量的reads，同时过滤去除多于5个N碱基的reads（-n 5）、总长度低于80个碱基的reads（-l 80），不启用滑窗统计功能（-M 0）。
构建参考基因组索引
给参考基因组构建索引
bwa index -a bwtsw Capsicum.annuum.L_Zunla-1_Release_2.0.fasta -p zunla
参考基因组比对与排序
1，比对
bwa mem -k 32 -M -t 4 -R '@RG\tID:S413\tPL:illumina\tLB:S413\tSM:S413' zunla S413.clean_R1.fastq.gz S413.clean_R2.fastq.gz | samtools view -S -b - > S413.bam
备注：zunla就是上面建立的参考基因组index，bwa mem输出的文件为sam文件，为文本文件，比较大，这里通过samtools工具把sam转为体积较小、读取较快的二进制bam文件。
2，比对完成之后对bam文件进行排序
samtools sort -@ 10 -m 80G -O bam -o S347.sorted.viewed.bam S413.bam
@是线程数，m参数是消耗的内存，我得是1t，所以随便设了
去除PCR扩增重复序列
由于构建文库的时候进行了pcr扩增，然后进行测序，这么测序出来的数据存在很多重复序列，在这里需要去除，然后建立索引，那么这个预处理工作就结束了。
遇到的问题最多的也是在这里。
按照之前的资料，对于PE测序数据一般用picard，效果好，SE可以用samtools或picard去重。
所以我首先用picard的MarkDuplicates进行了去重，遇到N个错误！
1，临时目录不够
错误提示：
Exception in thread "main" net.sf.picard.PicardException: Exception writing ReadEnds to file
Exception in thread "main" java.lang.OutOfMemoryError: GC overhead limit exceeded
这算是我节选的部分错误提示吧，看得出来，这是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
（2），为了一劳永逸的解决这个问题，我直接挂载了一个1t的分区给/tmp
有兴趣的可以借鉴我这个分区，给/tmp，/home和根目录/分别挂载1t，给数据目录/www挂载20t，还留几t机动调配。
2，java堆大小设置过低
错误提示：
Exception in thread "main" java.lang.OutOfMemoryError: Java heap space
面对“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”错误
INFO 2019-07-24 19:09:44 MarkDuplicates Read 704,000,000 records. Elapsed time: 01:21:51s. Time for last 1,000,000: 5s. Last read position: Chr12:227,643,683
INFO 2019-07-24 19:09:44 MarkDuplicates Tracking 5909309 as yet unmatched pairs. 32308 records in RAM.
[Wed Jul 24 19:12:35 CST 2019] picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 84.73 minutes.
Runtime.totalMemory()=54978936832
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 704691907, Read name A00174:10:HCCJKDSXX:4:2426:2564:33301, Insert size out of range
at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.java:455)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:849)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:834)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:802)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:574)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:553)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:543)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:257)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:305)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)
关于这个错误，从字面意思上来看是插入大小超过了范围，并且我发现基本都是在检测重复完成之后才出现，这个错误应当和程序读取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环境，所以上面的命令是我运行picard进行序列去重成功的命令。
如果不想使用conda的picard，那么我们也可以直接调用java包，在picard命令中已经包含了当前picard的jar包的地址，我们可以修改命令如下：
java -Xmx20g -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参数以后
我们可以看到大量的提示
Ignoring SAM validation error: ERROR: Record 691028838, Read name A00174:10:HCCJKDSXX:4:1515:25364:3208, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028841, Read name A00174:10:HCCJKDSXX:4:1352:30526:21104, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028854, Read name A00174:10:HCCJKDSXX:4:2563:12915:36949, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028860, Read name A00250:73:HFTMFDSXX:3:1456:13937:8563, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028895, Read name A00174:10:HCCJKDSXX:4:1435:10031:21746, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028896, Read name A00174:10:HCCJKDSXX:4:1572:29387:14653, Insert size out of range Ignoring SAM validation error: ERROR: Record 691028916, Read name A00174:10:HCCJKDSXX:4:1275:28519:35571, Insert size out of range Ignoring SAM validation error: ERROR: Record 691030765, Read name A00174:10:HCCJKDSXX:4:2545:1823:20588, Insert size out of range
其它去重软件
经过这一轮之后，整理发现，还有好几个很优秀的去重软件，
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)
这是由于我们的系统限制导致的，解决方法有两个。
（1），临时办法
首先
运行命令
ulimit -n
查看线程限制，默认一般为1024，我们修改为10240.
运行命令
ulimit -n 10240
（2），上面的方法重启后就失效了，如果想重启后生效，需要修改/etc/security/limits.conf文件，root权限。
在文件末尾添加
* soft nofile 10240 * hard nofile 10240
*代表所有该系统的用户
这个10240值不能大于/proc/sys/fs/nr_open文件中的数值，否则将无法正常登录系统。
我的nr_open值为1048576
2，samtools
samtools也有标记或者删除重复序列的功能，命令分别为
samtools rmdup
删除重复序列，但是软件还是推荐标记重复序列
samtools markdup -@ 8 S347.merged.sorted.viewed.bam S347.sorted.viewed.markdupba.bam
通过资料查询，发现这个工具对于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大小系统，并且最大参考序列长度为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 #去除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
