运行GATK时，报错：java.lang.IllegalArgumentException: samples cannot be empty

这个问题在GATK Forum有讨论到；

IllegalArgumentException: samples cannot be empty

Picard ValidateSamFile

ValidateSamFile下有两种mode:

VERBOSE : 检查到100个错误之后退出，并且输出错误到终端；

SUMMARYL: 输出结果是一个表格，展示errors和warnings的数目；

问题在bam文件 MISSING_READ_GROUP；请自动屏蔽WARNING结果；

bwa 比对使用参数-R

bwa -R 参数的使用

英文原版查看

测序时：

样本建一个库在同一条lane上完成测序产生reads sets可定义为一个Read group

样本建成分开独立的库测序得到的reads sets也就分属于不同的Read groups

存在于SAM/BAM /CRAM 文件中Read groups是由一系列标签组成；这些标签代表着测序中的一些技术特征；有了这些数据之后，大家就可以对bam文件进行重复序列标识和碱基质量重新矫正。

GATK要求输入的bam文件包含Read groups，如果没有就会报错。

例子：

GATk 要求read group的格式

  Read group是@RG开始。

ID = Read group identifier

  每一个Read group独有的ID；

  Illumina 测序数据中，read group IDs由flowcell ，lane name 和number组成。

  在矫正碱基质量时，read group IDs对区分技术批次效应是必须的；在这过程中，同一read group的reads假 定为有一样的技术误差。

PU = Platform Unit

  Platform Unit由三部分组成： {FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}

  {FLOWCELL_BARCODE} refers to the unique identifier for a particular flow cell

  The {LANE} indicates the lane of the flow cell

  The {SAMPLE_BARCODE} is a sample/library-specific identifier

  GATK 使用时，PU不是必须要求的；但是PU与ID同时存在时，PU优先级高于ID。

SM = Sample

  reads属于的样品名；SM要设定正确，因为GATK产生的VCF文件也使用这个名字。

PL = Platform/technology used to produce the read

  测序使用的平台： ILLUMINA, SOLID, LS454, HELICOS and PACBIO。

LB = DNA preparation library identifier

  对一个read group的reads进行重复序列标记时，需要使用LB来区分reads来自那条lane有时候，同一个库可能在不同的lane上完成测序为了加以区分，同一个或不同库只要是在不同的lane产生的reads都要单独给一个ID。

从read names中提取ID与PU

例子-Multi-sample and multiplexed example

三个样品：MOM, DAD, KID；

建库：每个样品分别建两个库，一个insert为200bp，一个insert为400bp；

上机：每个测序库使用Illumina HiSeq的 两条lane；

reads：来自 3 x 2 x 2 = 12条lane,可以产生12个bam文件，结果如下：

参考：

Read groups

Picard

从零开始完整学习全基因组测序（WGS）数据分析：第4节 构建WGS主流程

我们这里将用于流程构建的BWA就是其中最优秀的一个，它将BW(Burrows-Wheeler)压缩算法和后缀树相结合，能够让我们以较小的时间和空间代价，获得准确的序列比对结果。

别人的已安装文件打包传递后使用：

1、连接服务器

2、家目录下，Users，ls -all (或者 打开/etc里profile）

3、vim .bashrc

```

export PATH="your path:$PATH"

```

4、source .bashrc 

5、在全局中输入bwa测试是否可以使用

自己安装的软件路径都已经配置好了，所以不需要这样操作就可以直接使用。

BWA MEM比对模块是有一定适用范围的：它是专门为长read比对设计的，目的是为了解决，第三代测序技术这种能够产生长达几十kb甚至几Mbp的read情况。一般只有当read长度≥70bp的时候，才推荐使用，如果比这个要小，建议使用BWA ALN模块。

BWA ALN模块

1: 建立 Index

根据reference genome data(e.g. reference.fa) 建立 Index File

    bwa index -a bwtsw reference.fa

(bwa index 可以查看更多使用方法）

step 2: 寻找 SA coordinates

如果是pair-end 数据（leftRead.fastq和rightRead.fastq）两个文件分别处理

    bwa aln reference.fa leftRead.fastq >leftRead.sai

    bwa aln reference.fa rightRead.fastq >rightRead.sai

    bwa aln reference.fa singleRead.fastq >singleRead.sai

如果希望多线程运行，在其中加入 -t这个参数，另外-f这个参数可以指定结果输出文件，如:

    bwa aln -c -t 3 -f leftreads.sai reference.fa leftreads.fastq

step 3:转换SA coordinates输出为sam

如果是pair-end数据

    bwa sampe -f pair-end.sam reference.fa leftRead.sai rightRead.sai leftRead.fastq rightread.fastq

如果是single reads数据

    bwa samse -f single.sam reference.fa single.sai single.fastq

https://www.jianshu.com/p/859c0345624c

bwa mem [options] <idxbase><in1.fq>[in2.fq]

bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz >sample_name.sam

-t，线程数，我们在这里使用4个线程；-R 接的是 Read Group的字符串信息，这是一个非常重要的信息 ，以@RG开头，它是用来将比对的read进行分组的。不同的组之间测序过程被认为是相互独立的，这个信息对于我们后续对比对数据进行错误率分析和Mark duplicate时非常重要。在Read Group中，有如下几个信息非常重要：

(1) ID，这是Read Group的分组ID，一般设置为测序的lane ID（不同lane之间的测序过程认为是独立的），下机数据中我们都能看到这个信息的，一般都是包含在fastq的文件名中；

(2) PL， 指的是所用的测序平台，这个信息不要随便写！ 特别是当我们需要使用GATK进行后续分析的时候，更是如此！这是一个很多新手都容易忽视的一个地方，在GATK中，PL只允许被设置为：ILLUMINA，SLX，SOLEXA，SOLID，454，LS454，COMPLETE，PACBIO，IONTORRENT，CAPILLARY，HELICOS或UNKNOWN这几个信息。基本上就是目前市场上存在着的测序平台，当然，如果实在不知道，那么必须设置为UNKNOWN，名字方面不区分大小写。如果你在分析的时候这里没设置正确，那么在后续使用GATK过程中可能会碰到类似如下的错误。

(3) SM，样本ID，同样非常重要，有时候我们测序的数据比较多的时候，那么可能会分成多个不同的lane分布测出来，这个时候SM名字就是可以用于区分这些样本。

(4) LB，测序文库的名字，这个重要性稍微低一些，主要也是为了协助区分不同的group而存在。文库名字一般可以在下机的fq文件名中找到，如果上面的lane ID足够用于区分的话，也可以不用设置LB；

除了以上这四个之外，还可以自定义添加其他的信息，不过如无特殊的需要，对于序列比对而言，这4个就足够了。这些信息设置好之后， 在RG字符串中要用制表符（\t）将它们分开 。

sam转换bam（sam特殊二进制）

bwa mem -t 4 -R'@RG\tID:foo_lane\tPL:illumina\tLB:library\tSM:sample_name'/path/to/human.fasta read_1.fq.gz read_2.fq.gz | samtools view -S -b - >sample_name.bam

time$bwamem -t 40 -M -Y -R"@RG\tID:$RGID\tPL:ILLUMINA\tSM:$sample"$reference\$outdir/cleanfq/${fq_file_name}.paired.1.fq.gz$outdir/cleanfq/${fq_file_name}.paired.2.fq.gz |$samtoolsview -Sb - >$outdir/bwa/${sample}.bam &&\echo"** BWA MEM done **"&&\

-M 将 shorter split hits 标记为次优，以兼容 Picard’s markDuplicates 软件。

-Y Use soft clipping CIGAR operation for supplementary alignments. By default, BWA-MEM uses soft clipping for the primary alignment and hard clipping for supplementary alignments.

链接：https://www.jianshu.com/p/d6cc153be28d

设置权限 chmod 777 ./*.sh

运行

以上，我们就完成了read比对的步骤。接下来是排序：

-@，用于设定排序时的线程数；-m，每个线程排序时最大的内存消耗；-O 指定输出为bam格式；-o 是输出文件的名字.

在排序完成之后我们就可以开始执行去除重复（准确来说是  去除PCR重复序列 ）的步骤了。

把参数REMOVE_DUPLICATES设置为ture，那么重复序列就被删除掉，不会在结果文件中留存。我比较建议使用第一种做法，只是标记出来，并留存这些序列，以便在你需要的时候还可以对其做分析。

这一步完成之后，我们需要为sample_name.sorted.markdup.bam创建索引文件，它的作用能够让我们可以随机访问这个文件中的任意位置，而且后面的“局部重比对”步骤也要求这个BAM文件一定要有索引，命令如下：

$ samtools index sample_name.sorted.markdup.bam

接下来是局部区域重比对，通常也叫Indel局部区域重比对。有时在进行这一步骤之前还有一个merge的操作，将同个样本的所有比对结果合并成唯一一个大的BAM文件【注】，merge的例子如下：

$ samtools merge <out.bam><in1.bam>[<in2.bam>... <inN.bam>]

【注意】之所以会有这种情况，是因为有些样本测得非常深，其测序结果需要经过多次测序（或者分布在多个不同的测序lane中）才全部获得，这个时候我们一般会先分别进行比对并去除重复序列后再使用samtools进行合并。

---