PAKAP

INTRODUCTION (from paper?)

We have established a Present/Absent Kmer Analysis Pipeline (PAKAP). By reducing each read to component k-mers and comparing the relative abundance of these sub-sequences, we overcome statistical limitations of whole read comparative analysis.

PAKAP consists of a series of scripts written in Perl, Python and Bash scripts and requires Jellyfish [Marcais 2011] as well as optionally SOAPaligner. The scripts are freely available for non-commercial use.

What does PAKAP depend on?

• Jellyfish for fast kmer counting
• Some non-standard Perl modules:
  • bioperl
    • Bio::SeqIO
    • Bio::SearchIO
  • Parallel::ForkManager
  • Statistics::Descriptive
  • Config::IniFiles
  • GD::Graph::linespoints (for the script identifyKmerSize)
• Optional: SOAPaligner

Download

• Latest Version 1.0:
  • INSERT LINK

How to install?

• Download the [link_here DiffKAP package].

How to prepare the pipeline

Create your project configuration file by using the example config file. Here it is:

[PIPELINE]
# This should be the path to the directory containing the pipeline scripts.
SCRIPTS_PATH = /scratch/pipeline
# OPTIONAL: Path to reference
REFERENCE_PATH = /scratch/temp/ref.fasta
# OPTIONAL: Path to gene annotation on reference
GFF_PATH = /scratch/temp/ref.gff3

[JELLYFISH]
# The path to the local installation of jellyfish (http://www.cbcb.umd.edu/software/jellyfish/)
PATH = /scratch/bin
# The minimum k-mer size will be checked when KMER_SIZE=AUTO. Default is 5
MIN_KMER_SIZE=5
# The maximum k-mer size will be checked when KMER_SIZE=AUTO, jellyfish limits to the maximum k-mer size to 31. Default is 22
MAX_KMER_SIZE=22

[IDENTIFYKMERSIZE]
# Number of CPU will be used. Be aware of that the total amount of memory will be shared among all CPUs.
NUM_OF_PROCESSOR=4
# treatment 1 ID which will be used for naming files 
T1_ID=CON
# treatment 2 ID which will be used for naming files 
T2_ID=STE
# The output directory of the initial kmer counting, for jellyfish
OUT_DIR=/scratch/somewhere

# directory storing treatment 1 data files. The data files can be fasta or fastq formats.
DATA_DIR_T1=/scratch/somewhere/whole/whole_2.150.1sd/
# directory storing treatment 2 data files. The data files can be fasta or fastq formats.
DATA_DIR_T2=/scratch/somewhere/truncated/truncated_2.150.1sd/

[ADVANCED]
# Important setting in jellyfish for tuning jellyfish performance. A larger hash size, more memory will be used 
#   but less sub-count files will be generated. Default is 10000000.
JELLYFISH_HASH_SIZE=10000000
# Size of the jellyfish hash table. Use a size large enough to contain all of the K-mers such that 80% * s > number of distinct K-mers. [Default: 10000000].
#jellyfish_hash_size = 16G
# Length of counter in the hash table. See jellyfish manual for explanation. [Default: 4]
JELLYFISH_COUNTER_BITS = 18
# Number of bases used for splitting the k-mer files into sub-files. A larger number reduces RAM usage. Maximum at the moment is 3. [Default: 2]
N_SPLIT_BASES = 3
# Minimum occurrence count for a k-mer to be considered as a candidate Presence / Absence K-mer. [Default: 4]
MIN_OCC = 4

The most important settings are:

• SCRIPTS_PATH - this is where you extracted the pipeline
• REFERENCE_PATH - optional: if you have a reference fasta. SOAPaligner will be run if this is specified.
• GFF_PATH - optional: The SOAP results will be compared against genes if this is specified.
• JELLYFISH / PATH - the directory in which Jellyfish is located
• IDENTIFYKMERSIZE / OUT_DIR - the dirctory in which temporary kmers are stored (might get huge)
• DATA_DIR_T1 and DATA_DIR_T2 - where the fasta files for kmer counting are located

How to run the pipeline

The pipeline takes four arguments: 's1' specifies reads from the first group (or treatment etc.), 's2' specifies reads from the second group, 'c' is the path to the configuration file, and 'o' specifies the path to the directory which will contain all results.

The command is:

   python pipeline.py 
        --s1 ./truncated/truncated_2.150.1sd/truncated_2.150.1sd.R1.fasta ./truncated/truncated_2.150.1sd/truncated_2.150.1sd.R2.fasta 
        --s2 ./whole/whole_2.150.1sd/whole_2.150.1sd.R1.fasta ./whole/whole_2.150.1sd/whole_2.150.1sd.R2.fasta 
        -c default.config  
        -o output_folder/

or, easier:

   python pipeline.py 
        --s1 ./truncated/truncated_2.150.1sd/truncated_2.150.1sd.R?.fasta
        --s2 ./whole/whole_2.150.1sd/whole_2.150.1sd.R?.fasta
        -c default.config  
        -o output_folder/

This will read the configuration from default.config and generate all output files in output_folder.

How to interpret the results?

• You can download the results of the sample data here.

Reference

• Marçais, G. and Kingsford, C. (2011) A fast, lock-free approach for efficient parallel counting of occurrences of k-mers, Bioinformatics, 27, 764-770.

Back to Main_Page