-
Notifications
You must be signed in to change notification settings - Fork 2
Barcoding tools
The Barcoding tools are a little different from the other tools in that the do not access the BAM file directly, but rather they use VCFs which are created from BAM files from other tools. Hence the process is as follows:
- For each sample, we run
samtools mpileup(from samtools v1.3) to create a VCF file containing the allele read depths at each SNP in a master list (the SNPs we want in the barcode) - We then run the Barcoding tools to use the read counts in the VCFs to produce genotypes and then concatenate them in the desired order to form barcodes.
We need a list of SNPs that are concatenated to form the barcode. Indeed the list has to be in two different forms:
-
<SNPLIST_NAME>.tab- a tab-separated table of barcoding positions. This has one header line, then one line per barcoding position, in the desired barcode order. The table must have at least two columns, named "Chr" and "Pos", containing the chromosome ID and the position of the barcode variation within the chromosome, respectively. Any additional columns are ignored. The following is an example:
Num Chr Pos Ref Nonref 1 Pf3D7_02_v3 376222 A G 2 Pf3D7_02_v3 470013 G A 3 Pf3D7_03_v3 656861 T G [...]
-
<SNPLIST_NAME>.bed- a BED file of barcoding positions. This is needed bysamtoolsand is encoded in the BED format: a space-separate table (no header line), one line per barcode variation, ordered by genomic position. Each line contains three fields: the chromosome ID, and the start and end positions of the barcode variation. For SNPs, the end position is the same as the position in the<SNPLIST_NAME>.tabfile (see above), while the start position is one position earlier. [NOTE: as soon as I get round to it, I'll make a tool to generate the BED file from the SNP list automatically]. The following is an example of BED file, corresponding to the SNP list file shown above:
Pf3D7_02_v3 376221 376222 Pf3D7_02_v3 470012 470013 Pf3D7_03_v3 656860 656861 [...]
We produce one VCF file per sample using samtools, but then the file is reheadered (to fix the sample name) using bcftools. The commands are as follows:
samtools mpileup -I -R -C50 -d 1000 -f <REF_FASTA_FILE> -l <BED_FILE> -v -t AD,DP \ -o <RAW_VCFFILE> <BAM_FILE> bcftools reheader -s <SAMPLE_NAMEFILE> -o <FINAL_VCFFILE> <RAW_VCFFILE>
with the following parameter:
- REF_FASTA_FILE The path to a FASTA file containing the reference sequence used for alignment.
- BED_FILE The path to the BED file (see above).
- BAM_FILE The path to the BAM file containing the sample alignment.
- RAW_VCFFILE The path to the VCF file that will be output by mpileup.
-
FINAL_VCFFILE The path to the final, reheadered VCF file, e.g.
<SAMPLE_NAME>.vcf.gz - SAMPLE_NAMEFILE The path to a text file containing a single line with the sample name, e.g. "PA0007-C", used for reheadering
The final VCF has to be placed in an output folder <OUT_DIR>/<SUB> where <SUB> is a string consisting of the first 4 letters of the sample name. For example, file PA0001-C.alleles.tab will be written to folder <OUT_DIR>/PA00. This simple file hashing scheme, which is particularly suited for the MalariaGEN naming scheme, prevents thousands of files accumulating in a single folder, which may cause indexing problems especially on Linux.
The following script can be used as a template for creating the VCF file.
BAMLIST_FILE= BED_FILE= SNPLIST_FILE= REF_FASTA_FILE= OUT_DIR=
JOB=$LSB_JOBINDEX
IN=
sed "$JOB q;d" $BAMLIST_FILEread -a LINE <<< "$IN" SAMPLE=${LINE[0]} BAM_FILE=${LINE[1]}SUBDIR=$OUT_DIR/${SAMPLE:0:4} mkdir -p $SUBDIR
RAW_VCFFILE=$SUBDIR/$SAMPLE.raw.vcf.gz FINAL_VCFFILE=$SUBDIR/$SAMPLE.vcf.gz SAMPLE_NAMEFILE=$SUBDIR/$SAMPLE.sampleList.txt echo $SAMPLE > $SAMPLE_NAMEFILE
echo Starting $SAMPLE on $BAM_FILE samtools mpileup -I -R -C50 -d 1000 -f $REF_FASTA_FILE -l $BED_FILE -v -t AD,DP -o $RAW_VCFFILE $BAM_FILE bcftools reheader -s $SAMPLE_NAMEFILE -o $FINAL_VCFFILE $RAW_VCFFILE rm $SAMPLE_NAMEFILE rm $RAW_VCFFILE