【3.3.3】bam文件header中的@RG
比对软件bwa、bowtie2等产生的bam文件用于gatk分析时经常报错，提示没有read group信息。
Read group信息通常用来表示这些reads来自于同一个芯片的同一个run（个人理解：测序时试剂分别流过不同的lane，因此同一lane中的测序条件是相同的，排除了系统偏差）。最简单的例子：一个单独样品构建了一个文库，在一个flowcell的一个lane上进行测序，那么，最后得到的reads就是属于同一个group。实际上，一个lane可以测30G的数据，一般我们数据量小的话根本用不满一个lane，而是一个lane跑多个样品。现在下机数据已经根据index将同一lane中的不同样品reads分别存放到不同的文件中，因此同一文件中的reads都是同一group的。
一、基本信息
查看bam文件是否含有read group信息行：
samtools view -H sample.bam | grep '^@RG'
read group的格式（以gatk官方文档示例为例）：
@RG ID:H0164.2 PL:illumina PU:H0164ALXX140820.2 LB:Solexa-272222 PI:0 DT:2014-08-20T00:00:00-0400 SM:NA12878 CN:BI
以@RG开头，各个组成部分：
-
ID = read group identifier 必须要有。一个bam文件中ID必须是独一无二的，一般用“flowcell编号.lane编号”表示（示例中H0164即为flowcell编号，2为lane编号）。因此，同一文件中的reads来源于同一个run，可以用于后续base quality score recalibration（gatk3中关键一步，gatk4好像已经整合进HaplotypeCaller）等分析。
-
PU = platform unit 非必须。由{FLOWCELL_BARCODE}.{LANE}.{SAMPLE_BARCODE}构成，比ID多了样本的信息，优先级高于ID。
-
SM = sample 必须且唯一。表示该read group所存储的reads所来源的样本名称。同时也被用于VCF文件中的样本那一列。如果是测序了混合样品，则用pool name来代替。
-
PL = platform 必须。测序平台，只能是：ILLUMINA, SOLID, LS454, HELICOS and PACBIO。
-
LB = DNA preparation library identifier 必须。当一个DNA文库在不同lane上测序时，这个参数对于下游分析中的MarkDuplicates（标注建库过程中造成的duplicates）至关重要。
二、如何获得read group各项信息
测序平台、样本名称、文库名称这些信息测序人员已经完全掌握，直接提供即可。ID信息则可以通过下机数据fastq文件中的read名称得到。以个人通过Miseq平台测序得到数据为例：
A01415:337:HJFFGDSX3:4:1102:22571:31438
其中，A01415为机器编号，337为run id，HJFFGDSX3为flowcell编号，4为lane编号，1102为tile编号，22571为桥式扩增后的簇（cluster）在tile上的横坐标，31438为cluster在tile的纵坐标（一直误以为一个tile就是一个cluster。但是从这里可以看出，一个tile上有众多cluster）。
所以直接用ID:HJFFGDSX3.4作为ID即可。
上文说到一个fastq文件（双端测序可能是两个）中通常只有一个文库的一个lane的reads，已经满足了ID唯一的要求了。如果不放心，可以查看：
samtools view filename.bam | cut -f1 | cut -d ":" -f 3,4 | sort | uniq
三、修改
比对后修改RG
sample_name=SRR23455317;
time singularity exec /data2/neoantigen/singularity/nextNEOpi_1.3.2_18734d43.sif gatk AddOrReplaceReadGroups INPUT=work/"$sample_name".sam OUTPUT=work/"$sample_name".named.sam RGID=1 RGLB=lib1 RGPL=illumina RGPU=unit1 RGSM=1
参考资料
