annotate README @ 1:4f6952e0af48 default tip

CREST - add crest.loc.sample

This document has the information on how to run CREST for structural
variation detection.

=============
Requirements:
=============
Before running CREST, you need to make sure that several pieces of software
and/or modules are installed on the system:
1. BLAT software suite, especially blat, gfClient, and gfServer. BLAT
can be obtained from these links:
    BLAT for academic use: http://www.soe.ucsc.edu/~kent
    BLAT commercial license: http://www.kentinformatics.com/
2. CAP3 assembly program, available here:
    CAP3 for academic use: http://seq.cs.iastate.edu/cap3.html
    CAP3 commercial license: Contact Robin Kolehmainen at Michigan Tech, 
    [e-mail] or [phone].
3: SAMtools library for accessing SAM/BAM files, available from SourceForge:
    SAMtools: http://sourceforge.net/projects/samtools/files/
4. BioPerl and Bio::DB::Sam modules.  They are usually available as
packages on most Linux distributions, but are also available at this link:
    BioPerl: http://www.bioperl.org/
	Bio::DB::Sam: http://search.cpan.org/~lds/Bio-SamTools/lib/Bio/DB/Sam.pm
Important: you must install SAMtools library before install Bio::DB::Sam.
5. ptrfinder is needed if you want to remove short tandem repeat mediated
SVs, the executable is included in the download package, put it on the path.

Note:
1. You can use your own programs in place of BLAT and CAP3, but you need to
implement the run method in SVExtTools.pm.
2. The pipeline uses gfServer to mimic a standard blat server, so you need
to setup your own gfServer.  Details on setting up the server can be found in
the BLAT package.  Using a query server can significantly increase the speed of
the pipeline.

Your BAM files must contain soft-clipping signatures at the breakpoints.  If
they do not, you will not get any results.  For more information see the
section "About Soft-Clipping" at the end of this document.

=====================
Running the pipeline:
=====================

Make sure that all the required perl modules are in @INC.  One simple way
is to put all .pm and .pl scripts in the same directory and run them from
this same directory. Also, the input bam file, must be sorted and indexed
before running the pipeline.

We are going to use two sample bam files (tumor.bam and germline.bam) to
illustate how to run the pipeline. The examples assume you want to find SV
in tumor.bam and you also have the matched germline sample bam file.

Important: indexing all bam files before running the pipeline is required.

1. Get soft-clipping positions.

The program extractSClip.pl will extract all soft-clipping positions first,
and identify those positions with a cluster of soft-clipped reads.  The
program requires only the BAM file and the reference genome's FASTA file
The following is an example to extract all positions:

	extractSClip.pl -i tumor.bam --ref_genome hg18.fa

Two files named tumor.bam.cover and tumor.bam.sclip.txt will be generated
for use in the next step.

Note: The program can use either paired or single-end sequencing data.  For
single-end data, use the --nopaired parameter.

For whole genome sequencing project, we highly suggest running the procedure
in parallel by dividing the genome into pieces.  One natural way is by
chromosome. The following is an example to extract all positons on chr4.

	extractSClip.pl -i tumor.bam --ref_genome hg18.fa -r 4 

Important: The genome file used in this pipeline must be the same as the one
used to map reads, so the chromosome names need to agree.  In this example,
the genome file and bam file all have the chromosome name as 4 instead of
chr4 you may encounter.

Two files named tumor.bam.4.cover and tumor.bam.4.sclip.txt will be generated
for use in the next step. So it's very easy to run this step in parallel and
combine the results together to form a final result.

The output files for this step have names with suffixes of *.cover and
*.sclip.txt.  The .cover file is a tab-delimited text file, with columns:
chr, position, strand, number of soft-clipped reads, and coverage at that
position.  The strand is just left-clipped or right-clipped to help identify
the SV orientation.  The .sclip.txt file has the detailed information for
all soft-clipped reads including sequence and quality values. This file is
also tab-delimited with the following columns: chr, posiiton, strand, read
name, sequence, and quality.

Example of part of a *.cover file:
4       125892327       +       1       28
4       125892458       +       1       27
4       125893225       +       1       28
4       125893227       +       5       29
4       125893365       -       1       26
4       125893979       -       1       16
10      66301086        -       1       33
10      66301858        +       4       14
10      66301865        -       8       21
10      66301871        -       1       22
10      66302136        +       1       51

Example of part of a *.sclip.txt file:
4       125892327       +       HWUSI-EAS1591_6113C:3:17:12332:19420#0  CC      CC
4       125892458       +       HWUSI-EAS1591_6113C:4:91:6281:9961#0    GACTAACCACCACGGTACATGTTTTCCTATGTAAAAAACCTGCACATTCTACACATGTATCCCAGAACTTAAAGTAAAACAC B@C@?:CC>CCBCCCCACBCDCCCCCC;<:<9CCCCC@CCCCCBCCCCCCCCCCCCCCCACCCCCCCCCCCCCCCCCCCCAC
4       125893225       +       HWI-EAS90_614M9:5:18:17924:10181#0      CCCTCCTGGGTTCAAGTGATTCTCCTGCCTCTACCTCCCGAGTAGCTGGGATTACAGGTGCCCACCACCATGCCTGGCTAA  #######@@7@:8@><16+6(B>AABCAA3AB@CC6CCCCCCCDCCCCCCBCCDCCCCCCCCCCCCDCCCCCCCCCCCCCC

If you run this step in parallel, you need to combine the outputs by
concatenating the files.  Tumor and germline files must be concatenated
separately, for example:
cat tumor.bam.*.cover > tumor.bam.cover
cat germline.bam.*.cover > germline.bam.cover

2.Remove germline events (optional) 

Running step 1 on both germline and tumor samples, you will get the
soft-clipping posiitons in both samples.  This step will remove any position
in the tumor sample that also appears in germline sample, so germline events
will be removed.  This step does not use any sequence information and could
remove true events.  By our observations, true events are rarely removed.
You can skip this step and the program will do germline clean up at later step
(see the -g parameter for CREST.pl).

The script for this step is countDiff.pl and it only requires two parameters
to specity the two output files from previous step.

	countDiff.pl -d tumor.bam.cover -g germline.bam.cover > soft_clip.dist.txt

A file named tumor.bam.cover.somatic.cover will be generated for next step.

The program will generate a file with suffix *.somatic.cover, and it will
be used for the next step.  The file has the same format as *.cover generated
in the previous step.

The standard output will show the coverage distribution.  For every read count
in the range 1-999, it will show the number of breakpoints supported by that
many soft-clipped reads.

3. Running the SV detection script.

This is the core step in the detection process. The program is CREST.pl.

The program needs quite a few parameters, but you can think about what you
will need.  Here is a partial list of required and common parameters:

	-f The input soft-clipped coverage file produced in step 1 or 2.
	-d The disease or tumor bam file
	-g The germline bam file.  If you want to identify somatic SVs only, you
	should provide this parameter.  If you also want to identify germline
	events, you can leave this parameter unspecified. When treat your
	germline file as disease without specify -g parameter, the program can
	be used to identify germline events, or SV polymorphism.
	--ref_genome The reference genome in fa format (used by bam file)
	-t The reference genome in 2bit format (used by gfClient), this file can
	be generated by using faToTwoBit program in BLAT program suit.  This
	file must be the same as the one you used to setup gfServer.
	--blatserver The name or IP address of blat server, you need to use your
	own one instead of using the public one at UCSC.
	--blatport The port number for the blat server.
	--nopaired Tell the program the reads are not paired.

If all of the required programs are on the path then you won't need to
specify them again, otherwise you need to specify the paths to the programs
using the corresponding parameters.  Please use CREST.pl --man to show the 
man page, which provides a detailed parameter list.

An example of running this step is:
	CREST.pl -f tumor.bam.cover -d tumor.bam -g germline.bam --ref_genome \
	hg18.fa -t hg18.2bit

There is also a -r parameter to specify the range to be searched and it's highly
recommended to run using -r as below:
	CREST.pl -f tumor.bam.cover -d tumor.bam -g germline.bam --ref_genome \
	hg18.fa -t /genome/hg18.2bit -r chr1
So it's very easy to run the program in parallel by spliting the genome into 
pieces.

The program will generate a *.predSV.txt file.  The filename will be the input
bam with .predSV.txt appended unless you specify the -p parameter.  Also the 
STDERR output has the full list of SVs, including rejected ones.  The output 
file *.predSV.txt has the following tab-delimited columns: left_chr, left_pos, 
left_strand, # of left soft-clipped reads, right_chr, right_pos, right_strand, 
# right soft-clipped reads, SV type, coverage at left_pos, coverage at 
right_pos, assembled length at left_pos, assembled length at right_pos,
average percent identity at left_pos, percent of non-unique mapping reads at
left_pos, average percent identity at right_pos, percent of non-unique mapping
reads at right_pos, start position of consensus mapping to genome,
starting chromosome of consensus mapping, position of the genomic mapping of
consensus starting position, end position of consensus mapping to genome,
ending chromsome of consnesus mapping, position of genomic mapping of
consensus ending posiiton, and consensus sequences.  For inversion(INV), the 
last 7 fields will be repeated to reflect the fact two different breakpoints 
are needed to identify an INV event.

Example of the tumor.predSV.txt file:
4	125893227	+	5	10	66301858	-	4	CTX	29	14	83	71	0.895173453996983	0.230769230769231	0.735384615384615	0.5	1	4	125893135	176	10	66301773	TTATGAATTTTGAAATATATATCATATTTTGAAATATATATCATATTCTAAATTATGAAAAGAGAATATGATTCTCTTTTCAGTAGCTGTCACCTCCTGGGTTCAAGTGATTCTCCTGCCTCTACCTCCCGAGTAGCTGGGATTACAGGTGCCCACCACCATGCCTGGCTAATTTT
5	7052198	-	0	10	66301865	+	8	CTX	0	22	0	81	0.761379310344828	0.482758620689655	0	0	1	5	7052278	164	10	66301947	AGCCATGGACCTTGTGGTGGGTTCTTAACAATGGTGAGTCCGGAGTTCTTAACGATGGTGAGTCCGTAGTTTGTTCCTTCAGGAGTGAGCCAAGATCATGCCACTGCACTCTAGCCTGGGCAACAGAGGAAGACTCCACCTCAAAAAAAAAAAGTGGGAAGAGG
10	66301858	+	4	4	125893225	-	1	CTX	15	28	71	81	0.735384615384615	0.5	0.889507154213037	0.243243243243243	1	10	66301777	153	4	125893154	TTAGCCAGGCATGGTGGTGGGCACCTGTAATCCCAGCTACTCGGGAGGTAGAGGCAGGAGAATCACTTGAACCCAGGAGGTGACAGCTACTGAAAAGAGAATCATATTCTCTTTTCATAATTTAGAATATGATATATATTTCAAAATATGATA

If there are no or very few results, there may be a lack of soft-clipping.  See
the section "About Soft-Clipping" at the end of this document.

4. Visulization of the detailed alignment at breakpoint (optional)

The bam2html.pl script builds an html view of the multiple alignment for
the breakpoint, so you can manually check the soft-clipping and other things.

	bam2html.pl -d diag.bam -g germline.bam -i diag.bam.predSV.txt \
		--ref_genome /genome/hg18 -o diag.bam.predSV.html

The output file is specified by -o option.

====================
About Soft-Clipping:
====================

CREST uses soft-clipping signatures to identify breakpoints.  Soft-clipping is
indicated by "S" elements in the CIGAR for SAM/BAM records.  Soft-clipping may
not occur, depending on the mapping algorithm and parameters and sometimes even
the library preparation.

With bwa sampe:
---------------

One mapping method that will soft-clip reads is bwa sampe (BWA for paired-end
reads).  When BWA successfully maps one read in a pair but is not able to map
the other, it will attempt a more permissive Smith-Waterman alignment of the
unmapped read in the neighborhood of the mapped mate.  If it is only able to
align part of the read, then it will soft-clip the portion on the end that it
could not align.  Often this occurs at the breakpoints of structural
variations.

In some cases when the insert sizes approach the read length, BWA will not
perform Smith-Waterman alignment.  Reads from inserts smaller than the read
length will contain primer and/or adapter and will often not map.  When the
insert size is close to the read length, this creates a skewed distribution
of inferred insert sizes which may cause BWA to not attempt Smith-Waterman
realignment.  This is indicated by the error message "weird pairing".  Often
in these cases there are also unusually low mapping rates.

One way to fix this problem is to remap unmapped reads bwasw.  To do this,
extract the unmapped reads as FASTQ files (this may be done with a combination
of samtools view -f 4 and Picard's SamToFastq).  Realign using bwa bwasw and
build a BAM file.  Then, re-run CREST on this new BAM file, and you may pick
up events that would have been missed otherwise.

With other aligners:
--------------------

Consult the documentation or mailing list(s) for your mapper to determine its
behavior with regard to soft-clipping.