I have analysed some pooled RAD-seq data, and there is a particular program I want to use but it requires (1) one population, and (2) only one chromosome, from a mpileup file.
This is a bit of a hassle given that with population genomics with reduced representation you analyse multiple populations and fragments together, so the mpileup file has multiple populations in it from many different loci. Extracting single chromosomes from an mpileup is possible because of the indexing in the file structure, but is it possible to pull out which reads belong to each population?
I have provided a small example for "chromosome" (in the case of RAD-seq, its a "fragment") E137_L96 at the first 3 bases. There are four populations in the sample.
Does anyone know if this is possible? I know technically .bam files have each population's mapping information, but that is not accepted in the program I want to use.
Alternatively, does anyone know any good, NGS software, that can estimate heterozygosity from pooled data? I know PoPoolation exists, but to my knowledge it doesn't estimate 'H'.
Did you use unique barcodes for each population? Or did you just use the same barcode and sequence together. If the prior, then they should be pulled apart during demuxing. If the latter, then there is no way to tell which molecules came from which sample.
@Michael: Each population was barcoded separately. I have no trouble analysing the data set for variant calling, say in VarScan (each population comes out uniquely in the analysis), but what I want to do is go into the raw file and extract specific populations from it (using a script).
Essentially, I want to reconstruct the mpileup file to exclude/include only populations I am interested in. This information is obviously encoded in the file, I just can find any info about how this is done.
ReadGroups will allow you to do what you want. You can find some info about it here:
https://www.biostars.org/p/43897/
You should have a ReadGroup assigned to each of your sample alignments. If you havent already added ReadGroups while aligning, you can do this with PicardTools AddOrReplaceReadGroups
samtools mpileup allows to filter/exclude for ReadGroups with the -G option.