We sent 600 samples of a plant species (no reference genome) for DNA-extraction, library prep and ddRadSeq (paired end, Illumina Hiseq). Firstly, 200 samples were sequenced as pilot, then the others 400 followed the same procedure. The raw reads were then:
Demultiplexed and adepter clipped
Filtered by restriction enzyme cut site
Merged and Clusterized using CD-HIT-EST to form the reference contigs
Aligned against the reference contigs (Bowtie2)
Variant discovery and SNP called using Freebayes.
Then, many filters were applied to remove less informative SNPs or samples:
- filter out snps with missingness (% of samples with missing data) above 10%
- filter out samples with missingness (% of snps with missing data) above 10%
- filter out snps with minor allele frequency (frequency of the less represented allele) below 5%
This resulted in roughly 500 samples kept and 1500 snps.
The results show a clear batch effect between the samples sequenced as pilot and those sequenced after. Particularly:
The PCA show a clear separation on the first axis (explaining 3% of the total variance) between pilot and rest of the samples.
The samples from the pilot show an overall higher level of heterozygosity (12 % in pilot against 9 % rest of samples). This increased observation of heterozygous loci is distributed across many loci.
I've spent a few weeks trying to identify a technical factor that could explain these differences. The most remarkable technical difference I found between the samples is:
I calculated the median read count for each sample, then compared it between pilot and non-pilot samples. The mean of this parameter is higher in the pilot samples and the variance is half in comparison to the other samples. This is shown in the attached figure (AB: pilot samples, CDEF: other samples).
My questions are:
- Can a differential sequencing depth cause an increase/decrease of heterozygous call in radseq? Is there a read-processing step that could overcome this problem?
- Can you see other possible causes of these batch effects?