Difference between revisions of "SkimGBS"

Revision as of 05:37, 9 September 2014

This page collects the scripts and a manual for our SkimGBS genotyping by sequencing pipeline.

What SkimGBS depends on

For the Perl-scripts, just BioPerl. For the Python-scripts, BioPython and at least Python 2.7. we haven't tested the pipeline much with Python 3 but it should work.

Download

• Latest Version 1.0:
  • PAKAP_1.7.zip

How to install

Download the package from above, then just extract the files to a directory of your choosing. Remember the directory, since you'll need that in the next step.

How to run SkimGBS

Data

As input data, you need BAM alignments for your two parental individuals, and one BAM alignment for each individual of the DH or RIL population.

As a general overview: First align all reads, then call SNPs for the parental individuals using SGSautoSNP, then call genotypes for the population using snp_genotyping_all.pl, then filter SNPs using removeMonomorphicSNPs.py, then impute using imputeFlapjackSNPs.py, then visualize using Flapjack: http://ics.hutton.ac.uk/flapjack/

Alignments

You can filter reads based on quality first, however, since SGSautoSNP doesn't allow heterozygous SNPs, receiving a false-positive SNP from two reads that have exactly the same sequencing error seems unlikely.

Important: SGSautoSNP determines which read comes from which cultivar by checking the read name and expects either to see a P (for paired) or S (for single) followed by the number of the library. In practise, it's enough to use P1- for both cultivars. Example: We used Tapidor and Ningyou, so we added 'TP1-' and 'NP1-' to both read collections. You can do that either before you run the alignments, or after the alignments, here it is after the alignments using SOAPaligner:

   sed 's/^/TP1-/' all.soap > all_renamed.soap

We use SOAPaligner with the option -r 0 to exclude homeologous SNPs from reads aligning in several places, alternatively, filter BWA or Bowtie for quality scores, something above 20 or 30 should work.

   samtools view -q 30 mappings.bam > mappings_filtered.bam

It is best to use only reads that align in pairs in order to increase the accuracy of called SNPs, so depending on your aligner you will have to filter unpaired reads. In SOAPaligner, just use the paired output files.

We merged all files for the two parents:

  samtools merge TapidorAndNingyou.bam Tapidor.bam Ningyou.bam

Remove clones using Picard's MarkDuplicates:

  java -jar `which MarkDuplicates.jar` INPUT=TapidorAndNingyou.bam OUTPUT=TapidorAndNingyou_cleaned.bam METRICS_FILE=stats.txt REMOVE_DUPLICATES=true

This overfilters a bit, since by chance you will find some reads align from different libraries in exactly the same position, so they may be treated as clones. Alternatively, you can run MarkDuplicates for each library itself and merge the results.

To make sure everything worked fine, use Picard again:

  java -jar `which ValidateSamFile.jar` INPUT=your_output_file

SGSautoSNP can either run one process per reference chromosome using multiple threads, or you can run many SGSautoSNP runs, one per chromosome. I prefer the latter, so I split the merged parental reads by chromosomes using bamtools. Both give you exactly the same results, running one SGSautoSNP per chromosome is easier to do when you have job arrays as with PBS.

  bamtools split -in TapidorAndNingyou_cleaned.bam -reference

This will result in one BAM file per chromosome which includes the data of both parental cultivars.

SNP-calling

SGSautoSNP can produce SNPs for FASTA files containing one or several references. To make sure that it offsets correctly, you have to supply a gff3 file of contigs, even if it's only one chromomosome.

FAQ

To come.

Reference

Back to Main_Page