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)
• Python 2.7, currently the pipeline won't run under Python 3
  • BioPython
• Optional: SOAPaligner

Download

• Latest Version 1.0:
  • INSERT LINK

How to install?

• Download the [link_here DiffKAP package].

How to prepare the pipeline

IMPORTANT: Currently, only single-line FASTA files are accepted by the pipeline, in other words, one line per identifier, one line per sequence. If you have FASTQ files, the included script FixFASTQ.py can transform them into the needed format. Usage:

   python FixFASTQ.py your_reads.fastq single_line_reads.fasta

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 [link_here here].

In the given out_dir, between two and four folders will be created.

1. kmers - these are the k-mers present in s1 only, and in s2 only. K-mers present in s1 only start with s1.

2. PARs - this folder contains the reads present only in one of the two conditions:

2.1 s1_0_PARs_A.fasta - the first partner reads of the pairs only present in s1

2.2 s1_0_PARs_B.fasta - the second partner reads of the pairs only present in s1

2.3 s1_0_PARs_single.fasta - the pairs, together with unpaired reads only present in s1

2.4 the same for s2.

3. SOAPs - optional, only when reference is specified: The SOAP results for aligning the paired and single reads from s1 and s2. Can be converted to SAM/BAM using SOAPaligner's soap2sam.pl and samtools view.

4. Results - optional, only when gff file and reference are specified: A table detailing how many genes there are in the given gff-file, how many reads were produced by the pipeline, how many reads could be aligned, how many reads hit inside the genes in the gff3 file, how many reads hit only with their start or their end inside genes, how many reads hit multiple times etc. Example:

Name	Number_of_genes	Aligned_reads	Complete_hit	Only_start	Only_end	Number_alignments_without_hit_to_gene	Duplicate_hits
chr1.232_genes_vs_truncated_10.100.1sd_s2_0	232	2670	2655	8	7	0	0

Here, the large majority of produced reads aligns completely inside the given 232 genes (2655 out of 2670 reads). 8 reads align only with their start and then stretch outside the gene, 7 reads start before the gene and then end inside the gene.

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