opensource/galaxy
3
0
Fork 0
mirror of https://github.com/galaxyproject/galaxy.git synced 2026-10-29 17:42:12 +08:00
Code Issues Packages Projects Releases Wiki Activity
Files
dffac945d2ff2ca5c29082bcb2bf4df4d5764955
galaxy/tools/sr_mapping/lastz_paired_reads_wrapper.py
T

848 lines
43 KiB
Python
Raw Blame History

#!/usr/bin/env python
"""
Runs Lastz paired read alignment process
Written for Lastz v. 1.02.00.
# Author(s): based on various scripts written by Bob Harris (rsharris@bx.psu.edu),
# then tweaked to this form by Greg Von Kuster (greg@bx.psu.edu)
This tool takes the following input:
a. A collection of 454 paired end reads ( a fasta file )
b. A linker sequence ( a very small fasta file )
c. A reference genome ( nob, 2bit or fasta )
and uses the following process:
1. Split reads into mates: the input to this step is the read file XXX.fasta, and the output is three
files; XXX.short.fasta, XXX.long.fasta and XXX.mapping. The mapping file records the information necessary
to convert mate coordinates back into the original read, which is needed later in the process.
2. Align short mates to the reference: this runs lastz against every chromosome. The input is XXX.short.fasta
and the reference genome, and the output is a SAM file, XXX.short.sam.
3. Align long mates to the reference: this runs lastz against every chromosome. The input is XXX.long.fasta
and the reference genome, and the output is a SAM file, XXX.long.sam.
4. Combine, and convert mate coordinates back to read coordinates. The input is XXX.mapping, XXX.short.sam and
XXX.long.sam, and the output is XXX.sam.
usage: lastz_paired_reads_wrapper.py [options]
--ref_name: The reference name to change all output matches to
--ref_source: The reference is cached or from the history
--source_select: Use pre-set or cached reference file
--input1: The name of the reference file if using history or reference base name if using cached
--input2: The reads file to align
--input3: The sequencing linker file
--input4: The base quality score 454 file
--ref_sequences: The number of sequences in the reference file if using one from history
--output: The name of the output file
--lastz_seqs_file_dir: Directory of local lastz_seqs.loc file
"""
import optparse, os, subprocess, shutil, sys, tempfile, time
from string import maketrans
from galaxy import eggs
import pkg_resources
pkg_resources.require( 'bx-python' )
from bx.seq.twobit import *
from bx.seq.fasta import FastaReader
from galaxy.util.bunch import Bunch
from galaxy.util import string_as_bool
# Column indexes for SAM required fields
SAM_QNAME_COLUMN = 0
SAM_FLAG_COLUMN = 1
SAM_RNAME_COLUMN = 2
SAM_POS_COLUMN = 3
SAM_MAPQ_COLUMN = 4
SAM_CIGAR_COLUMN = 5
SAM_MRNM_COLUMN = 6
SAM_MPOS_COLUMN = 7
SAM_ISIZE_COLUMN = 8
SAM_SEQ_COLUMN = 9
SAM_QUAL_COLUMN = 10
SAM_MIN_COLUMNS = 11
# SAM bit-encoded flags
BAM_FPAIRED = 1 # the read is paired in sequencing, no matter whether it is mapped in a pair
BAM_FPROPER_PAIR = 2 # the read is mapped in a proper pair
BAM_FUNMAP = 4 # the read itself is unmapped; conflictive with BAM_FPROPER_PAIR
BAM_FMUNMAP = 8 # the mate is unmapped
BAM_FREVERSE = 16 # the read is mapped to the reverse strand
BAM_FMREVERSE = 32 # the mate is mapped to the reverse strand
BAM_FREAD1 = 64 # this is read1
BAM_FREAD2 = 128 # this is read2
BAM_FSECONDARY = 256 # not primary alignment
BAM_FQCFAIL = 512 # QC failure
BAM_FDUP = 1024 # optical or PCR duplicate
# Keep track of all created temporary files so they can be deleted
global tmp_file_names
tmp_file_names = []
# The values in the skipped_lines dict are tuples consisting of:
# - the number of skipped lines for that error
# If not a sequence error:
# - the 1st line number on which the error was found
# - the text of the 1st line on which the error was found
# If a sequence error:
# - The number of the sequence in the file
# - the sequence name on which the error occurred
# We may need to improve dealing with file position and text as
# much of it comes from temporary files that are created from the
# inputs, and not the inputs themselves, so this could be confusing
# to the user.
global skipped_lines
skipped_lines = dict( bad_interval=( 0, 0, '' ),
inconsistent_read_lengths=( 0, 0, '' ),
inconsistent_reads=( 0, 0, '' ),
inconsistent_sizes=( 0, 0, '' ),
missing_mate=( 0, 0, '' ),
missing_quals=( 0, 0, '' ),
missing_seq=( 0, 0, '' ),
multiple_seqs=( 0, 0, '' ),
no_header=( 0, 0, '' ),
num_fields=( 0, 0, '' ),
reads_paired=( 0, 0, '' ),
sam_flag=( 0, 0, '' ),
sam_headers=( 0, 0, '' ),
sam_min_columns=( 0, 0, '' ),
two_mate_names=( 0, 0, '' ),
wrong_seq_len=( 0, 0, '' ) )
global total_skipped_lines
total_skipped_lines = 0
def stop_err( msg ):
sys.stderr.write( "%s" % msg )
sys.exit()
def skip_line( error_key, position, text ):
if not skipped_lines[ error_key ][2]:
skipped_lines[ error_key ][1] = position
skipped_lines[ error_key ][2] = text
skipped_lines[ error_key ][0] += 1
total_skipped_lines += 1
def get_tmp_file_name( dir=None, suffix=None ):
"""
Return a unique temporary file name that can be managed. The
file must be manually removed after use.
"""
if dir and suffix:
tmp_fd, tmp_name = tempfile.mkstemp( dir=dir, suffix=suffix )
elif dir:
tmp_fd, tmp_name = tempfile.mkstemp( dir=dir )
elif suffix:
tmp_fd, tmp_name = tempfile.mkstemp( suffix=suffix )
os.close( tmp_fd )
tmp_file_names.append( tmp_name )
return tmp_name
def run_command( command ):
proc = subprocess.Popen( args=command, shell=True, stderr=subprocess.PIPE, )
proc.wait()
stderr = proc.stderr.read()
proc.wait()
if stderr:
stop_err( stderr )
def split_paired_reads( input2, combined_linker_file_name ):
"""
Given a fasta file of allegedly paired end reads ( input2 ), and a list of intervals
showing where the linker is on each read ( combined_linker_file_name ), split the reads into left and right
halves.
The input intervals look like this. Note that they may include multiple intervals for the same read
( which should overlap ), and we use the union of them as the linker interval. Non-overlaps are
reported to the user, and those reads are not processed. Starts are origin zero.
#name strand start len size
FG3OYDA05FTEES + 219 42 283
FG3OYDA05FVOLL + 263 41 416
FG3OYDA05FFL7J + 81 42 421
FG3OYDA05FOQWE + 55 42 332
FG3OYDA05FV4DW + 297 42 388
FG3OYDA05FWAQV + 325 42 419
FG3OYDA05FVLGA + 90 42 367
FG3OYDA05FWJ71 + 58 42 276
The output gives each half-sequence on a separate line, like this. This allows easy sorting of the
sequences by length, after the fact.
219 FG3OYDA05FTEES_L TTTAGTTACACTTAACTCACTTCCATCCTCTAAATACGTGATTACCTTTC...
22 FG3OYDA05FTEES_R CCTTCCTTAAGTCCTAAAACTG
"""
# Bob says these should be hard-coded.
seq_len_lower_threshold = 17
short_mate_cutoff = 50
# We need to pass the name of this file back to the caller.
tmp_mates_file_name = get_tmp_file_name( suffix='mates.txt' )
mates_file = file( tmp_mates_file_name, "w+b" )
# Read the linker intervals
combined_linker_file = file( combined_linker_file_name, "rb" )
read_to_linker_dict = {}
i = 0
for i, line in enumerate( combined_linker_file ):
line = line.strip()
if line.startswith( "#" ):
continue
if line.find( '#' ) >= 0:
line = line.split( "#", 1 )[0].rstrip()
fields = line.split()
if len( fields ) != 4:
skip_line( 'num_fields', i+1, line )
continue
name, start, length, size = fields
start = int( start )
length = int( length )
size = int( size )
end = start + length
if end > size:
skip_line[ 'bad_interval' ] += 1
continue
if name not in read_to_linker_dict:
read_to_linker_dict[ name ] = ( start, end, size )
continue
if read_to_linker_dict[ name ] == None:
# Read previously marked as non-overlapping intervals, so skip this sequence - see below
continue
( s, e, sz ) = read_to_linker_dict[ name ]
if sz != size:
skip_line( 'inconsistent_sizes', i+1, name )
continue
if s > end or e < start:
# Non-overlapping intervals, so skip this sequence
read_to_linker_dict[ name ] = None
continue
read_to_linker_dict[ name ] = ( min( s, start ), max( e, end ), size )
combined_linker_file.close()
# We need to pass the name of this file back to the caller.
tmp_mates_mapping_file_name = get_tmp_file_name( suffix='mates.mapping' )
mates_mapping_file = file( tmp_mates_mapping_file_name, 'w+b' )
# Process the sequences
seqs = 0
fasta_reader = FastaReader( file( input2, 'rb' ) )
while True:
seq = fasta_reader.next()
if not seq:
break
seqs += 1
if seq.name not in read_to_linker_dict:
if seq.length > seq_len_lower_threshold:
mates_file.write( "%-3d %s %s\n" % ( seq.length, seq.name, seq.text ) )
read_to_linker_dict[ seq.name ] = ""
continue
if read_to_linker_dict[ seq.name ] == "":
skip_line( 'multiple_seqs', seqs, seq.name )
continue
if read_to_linker_dict[ seq.name ] == None:
# Read previously marked as non-overlapping intervals, so skip this sequence - see above
continue
( start, end, size ) = read_to_linker_dict[ seq.name ]
if seq.length != size:
skip_line( 'wrong_seq_len', seqs, seq.name )
continue
left = seq.text[ :start ]
right = seq.text[ end: ]
left_is_small = len( left ) <= seq_len_lower_threshold
right_is_small = len( right ) <= seq_len_lower_threshold
if left_is_small and right_is_small:
continue
if not left_is_small:
mates_file.write( "%-3d %s %s\n" % ( len( left ), seq.name + "_L", left ) )
mates_mapping_file.write( "%s %s %s %s\n" % ( seq.name + "_L", seq.name, 0, size - start ) )
if not right_is_small:
mates_file.write( "%-3d %s %s\n" % ( len( right ), seq.name + "_R", right ) )
mates_mapping_file.write( "%s %s %s %s\n" % ( seq.name + "_R", seq.name, end, 0 ) )
read_to_linker_dict[ seq.name ] = ""
combined_linker_file.close()
mates_file.close()
mates_mapping_file.close()
# Create temporary files for short and long mates
tmp_mates_short_file_name = get_tmp_file_name( suffix='mates.short' )
tmp_mates_long_file_name = get_tmp_file_name( suffix='mates.long' )
tmp_mates_short = open( tmp_mates_short_file_name, 'w+b' )
tmp_mates_long = open( tmp_mates_long_file_name, 'w+b' )
i = 0
for i, line in enumerate( file( tmp_mates_file_name, 'rb' ) ):
fields = line.split()
seq_len = int( fields[0] )
seq_name = fields[1]
seq_text = fields[2]
if seq_len <= short_mate_cutoff:
tmp_mates_short.write( ">%s\n%s\n" % ( seq_name, seq_text ) )
else:
tmp_mates_long.write( ">%s\n%s\n" % ( seq_name, seq_text ) )
tmp_mates_short.close()
tmp_mates_long.close()
return tmp_mates_mapping_file_name, tmp_mates_file_name, tmp_mates_short_file_name, tmp_mates_long_file_name
def align_mates( input1, ref_source, ref_name, ref_sequences, tmp_mates_short_file_name, tmp_mates_long_file_name ):
tmp_align_file_names = []
if ref_source == 'history':
# Reference is a fasta dataset from the history
# Create temporary files to contain the output from lastz executions
tmp_short_file_name = get_tmp_file_name( suffix='short_out' )
tmp_align_file_names.append( tmp_short_file_name )
tmp_long_file_name = get_tmp_file_name( suffix='long_out' )
tmp_align_file_names.append( tmp_long_file_name )
seqs = 0
fasta_reader = FastaReader( open( input1 ) )
while True:
# Read the next sequence from the reference dataset. Note that if the reference contains
# a small number of chromosomes this loop is ok, but in many cases the genome has a bunch
# of small straggler scaffolds and contigs and it is a computational waste to do each one
# of these in its own run. There is an I/O down side to running by subsets (even if they are
# one sequence per subset), compared to splitting the reference into sizes of 250 mb. With
# the subset action, lastz still has to read and parse the entire file for every run (this
# is true for fasta, but for .2bit files it can access each sequence directly within the file,
# so the overhead is minimal).
"""
:> output_file (this creates the output file, empty)
while there are more sequences to align
find the next sequences that add up to 250M, put their names in farf.names
lastz ${refFile}[subset=farf.names][multi][unmask] ${matesPath}/${matesFile} ...
>> output_file
"""
seq = fasta_reader.next()
if not seq:
break
seqs += 1
# Create a temporary file to contain the current sequence as input to lastz.
# We're doing this a bit differently here since we could be generating a huge
# number of temporary files.
tmp_in_fd, tmp_in_file_name = tempfile.mkstemp( suffix='seq_%d_in' % seqs )
tmp_in_file = os.fdopen( tmp_in_fd, 'w+b' )
tmp_in_file.write( '>%s\n%s\n' % ( seq.name, seq.text ) )
tmp_in_file.close()
# Align short mates
command = 'lastz %s[unmask]%s %s ' % ( tmp_in_file_name, ref_name, tmp_mates_short_file_name )
command += 'Z=1 --seed=1111111011111 --notrans --maxwordcount=90% --match=1,3 O=1 E=3 X=15 K=10 Y=12 L=18 --ambiguousn --noytrim --identity=95 --coverage=80 --continuity=95 --format=softsam- '
command += '>> %s' % tmp_short_file_name
run_command( command )
# Align long mates
command = 'lastz %s[unmask]%s %s ' % ( tmp_in_file_name, ref_name, tmp_mates_long_file_name )
command += 'Z=15 W=13 --notrans --exact=18 --maxwordcount=90% --match=1,3 O=1 E=3 Y=10 L=18 --ambiguousn --noytrim --identity=95 --coverage=90 --continuity=95 --format=softsam- '
command += '>> %s' % tmp_long_file_name
run_command( command )
# Remove the temporary file that contains the current sequence
os.remove( tmp_in_file_name )
else:
# Reference is a locally cached 2bit file, split lastz calls across number of chroms in 2bit file
tbf = TwoBitFile( open( input1, 'rb' ) )
for chrom in tbf.keys():
# Align short mates
tmp_short_file_name = get_tmp_file_name( suffix='short_vs_%s' % chrom )
tmp_align_file_names.append( tmp_short_file_name )
command = 'lastz %s/%s[unmask]%s %s ' % ( input1, chrom, ref_name, tmp_mates_short_file_name )
command += 'Z=1 --seed=1111111011111 --notrans --maxwordcount=90% --match=1,3 O=1 E=3 X=15 K=10 Y=12 L=18 --ambiguousn --noytrim --identity=95 --coverage=80 --continuity=95 --format=softsam- '
command += '> %s' % tmp_short_file_name
run_command( command )
# Align long mates
tmp_long_file_name = get_tmp_file_name( suffix='long_vs_%s' % chrom )
tmp_align_file_names.append( tmp_long_file_name )
command = 'lastz %s/%s[unmask]%s %s ' % ( input1, chrom, ref_name, tmp_mates_long_file_name )
command += 'Z=15 W=13 --notrans --exact=18 --maxwordcount=90% --match=1,3 O=1 E=3 Y=10 L=18 --ambiguousn --noytrim --identity=95 --coverage=90 --continuity=95 --format=softsam- '
command += '> %s' % tmp_long_file_name
run_command( command )
return tmp_align_file_names
def paired_mate_unmapper( input2, input4, tmp_mates_mapping_file_name, tmp_align_file_name_list, output ):
"""
Given a SAM file corresponding to alignments of *subsegments* of paired 'reads' to a reference sequence,
convert the positions on the subsegments to positions on the reads. Also (optionally) add quality values.
The input file is in SAM format, as shown below. Each line represents the alignment of a part of a read
to a reference sequence. Read pairs are indicated by suffixes in their names. Normally, the suffixes _L
and _R indicate the left and right mates of reads (this can be overridden with the --left and --right
options). Reads that were not mates have no suffix.
(SAM header lines omitted)
F2YP0BU02G7LK5_R 16 chr21 15557360 255 40M * 0 0 ATTTTATTCTCTTTGAAGCAATTGTGAATGGGAGTTTACT *
F2YP0BU02HXV58_L 16 chr21 15952091 255 40M6S * 0 0 GCAAATTGTGCTGCTTTAAACATGCGTGTGCAAGTATCTTtttcat *
F2YP0BU02HREML_R 0 chr21 16386077 255 33M5S * 0 0 CCAAAGTTCTGGGATTACAGGCGTGAGCCATCGcgccc *
F2YP0BU02IOF1F_L 0 chr21 17567321 255 7S28M * 0 0 taaagagAAGAATTCTCAACCCAGAATTTCATATC *
F2YP0BU02IKX84_R 16 chr21 18491628 255 22M1D18M9S * 0 0 GTCTCTACCAAAAAATACAAAAATTAGCCGGGCGTGGTGGcatgtctgt *
F2YP0BU02GW5VA_L 16 chr21 20255344 255 6S32M * 0 0 caagaaCAAACACATTCAAAAGCTAGTAGAAGGCAAGA *
F2YP0BU02JIMJ4_R 0 chr21 22383051 255 19M * 0 0 CCCTTTATCATTTTTTATT *
F2YP0BU02IXZGF_L 16 chr21 23094798 255 13M1I18M * 0 0 GCAAGCTCCACTTCCCGGGTTCACGCCATTCT *
F2YP0BU02IODR5_L 0 chr21 30935325 255 37M * 0 0 GAAATAAAGGGTATTCAATTAGGAAAAGAGGAAGTCA *
F2YP0BU02IMZBL_L 16 chr21 31603486 255 28M1D1M * 0 0 ATACAAAAATTAGCCGGGCACAGTGGCAG *
F2YP0BU02JA9PR_L 16 chr21 31677159 255 23M * 0 0 CACACCTGTAACCCCAGCACTTT *
F2YP0BU02HKC61_R 0 chr21 31678718 255 40M * 0 0 CACTGCACTCCAGCCTGGGTGACAAAGCAAGACTCTGTCT *
F2YP0BU02HKC61_R 0 chr21 31678718 255 40M * 0 0 CACTGCACTCCAGCCTGGGTGACAAAGCAAGACTCTGTCT *
F2YP0BU02HVA88 16 chr21 31703558 255 1M1D35M8S * 0 0 TGGGATTACAGGCGTGAGCTACCACACCCAGCCAGAgttcaaat *
F2YP0BU02JDCF1_L 0 chr21 31816600 255 38M * 0 0 AGGAGAATCGCTTGAACCCAGGAGGCAGAGGTTGCGGT *
F2YP0BU02GZ1GO_R 0 chr21 33360122 255 6S38M * 0 0 cctagaCTTCACACACACACACACACACACACACACACACACAC *
F2YP0BU02FX387_L 16 chr22 14786201 255 26M * 0 0 TGGATGAAGCTGGAAACCATCATTCT *
F2YP0BU02IF2NE_R 0 chr22 16960842 255 40M10S * 0 0 TGGCATGCACCTGTAGTCTCAGCTACTTGGGAGGCTGAGGtgggaggatc *
F2YP0BU02F4TVA 0 chr22 19200522 255 49M * 0 0 CCTGGGAGGCGGAGGTTGCAGTGAGCCGAGATCACGCCATTGCACTCCA *
F2YP0BU02HKC61_R 16 chr22 29516998 255 8S32M * 0 0 agacagagTCTTGCTTTGTCACCCAGGCTGGAGTGCAGTG *
F2YP0BU02FS4EM_R 0 chr22 30159364 255 29M * 0 0 CTCCTGCCTCAGCCTCCCGAGTAGTTGGG *
F2YP0BU02G197P_L 0 chr22 32044496 255 40M10S * 0 0 TTGTTGGACATTTGGGTTGGTTCCAAGTCTTTGCTATTGTgaataatgcc *
F2YP0BU02FIING 16 chr22 45959944 255 3M1I11M1I26M * 0 0 AGCTATGGTACTGGCTATGAAAGCAGACACATAGACCAATGG *
F2YP0BU02GUB9L_L 16 chr22 49198404 255 16M1I20M * 0 0 CACCACGCTCGGCTAATTTTTGTATTTTTAGTAGAGA *
The user must provide a mapping file (which might better be called an unmapping file). This file is usually
created by split_paired_reads, and tells us how to map the subsegments back to original coordinates in a single
read (this means the left and right mates were part of a single read). The mapping file contains four columns.
The first two give the mates's name (including the suffix) and the read name. The last two columns describe how
much of the full original sequence is missing from the mate. For example, in the read below, the left mate is
missing 63 on the right (42 for the linker and 21 for the right half). The right mate is missing 339 on the left.
left half: TTTCAACATATGCAAATCAATAAATGTAATCCAGCATATAAACAGAACCA
AAGACAAAAACCACATGATTATCTCAATAGATGCAGAAAAGGCCTTCGGC
AAAATTCAACAAAACTCCATGCTAAAACTCTCAATAAGGTATTGATGGGA
CATGCCGCATAATAATAAGACATATCTATGACAAACCCACAGCCAATATC
ATGCTGAATGCACAAAAATTGGAAGCATTCCCTTTGAAAACTGGCACAAG
ACTGGGATGCCCTCTCTCACAACTCCTATTCAACATAGTGTTGGAAG
linker: CGTAATAACTTCGTATAGCATACATTATACGAAGTCATACGA
right half: CTCCTGCCTCAGCCTCCCGAG
mate_name read_name offset_to_start offset_from_end
F2YP0BU02FS4EM_L F2YP0BU02FS4EM 0 71
F2YP0BU02FS4EM_R F2YP0BU02FS4EM 339 0
The user can also specify a quality scores file, which should look something like this. Quality values are presumed
to be PHRED scores, written in space-delimited decimal.
>F2YP0BU02FS4EM
38 38 38 40 40 40 40 40 40 40 40 40 40 39 39 39 40 40 40 40 38 21 21 21 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 39 39 39 40 40 40 40 40 40 40 33
32 32 40 40 40 21 21 18 18 21 34 34 31 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 32 32 32 32 40 40 40 40 40 40 40 34 34 35
31 31 28 28 33 33 33 36 36 36 17 17 17 19 26 36 36 36 40 40 40 40 40 33 34
34 34 39 39 39 40 40 40 40 40 33 33 34 34 40 40 40 40 40 40 40 39 39 39 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 39 39 39 39 39 39 40 40 40 39 39 39 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 40 26 26 26 26 26 40 40 38 38 37 35 33
36 40 19 17 17 17 17 19 19 23 30 20 20 20 23 35 40 36 36 36 36 36 36 36 36
39 40 34 20 27 27 35 39 40 37 40 40 40 40 40 40 40 40 40 40 34 34 35 39 40
40 40 40 40 40 40 39 39 39 40 40 40 40 36 36 32 32 28 28 29 30 36 40 30 26
26 26 34 39 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 39 39 39
40 39 35 34 34 40 40 40 40 30 30 30 35 40 40 40 40 40 39 39 36 40 40 40 40
39 39 39 39 30 30 28 35 35 39 40 40 40 40 40 35 35 35
>F2YP0BU02G197P
40 40 40 40 40 40 40 40 40 40 39 39 39 39 39 39 40 40 40 40 40 40 40 40 40
40 40 40 40 26 26 26 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 34 34 34 40 40
40 40 40 40 40 40 39 39 39 40 40 40 40 40 40 40 40 40 40 39 39 39 40 40 40
40 40 40 40 40 40 40 34 34 34 34 40 40 40 40 34 34 34 34 40 40 40 40 40 40
40 40 40 40 40 39 39 39 34 34 34 34 40 40 40 40 39 39 25 25 26 39 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
33 33 33 33 40 35 21 21 21 30 38 40 40 40 40 40 40 40 40 35 35 30 30 30 40
40 40 39 39 39 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40 40
40 40 40 40 40 40 40 40 40 40 40 40 39 39 39 40 40 40 40 40 40 40 40 40 40
40 40 40 39 39 39 40 40
>F2YP0BU02FIING
32 32 32 25 25 25 25 24 25 30 31 30 27 27 27 28 28 21 19 19 13 13 13 14 19
19 17 19 16 16 25 28 22 21 17 17 18 25 24 25 25 25
The output file is also SAM:
(SAM header lines omitted)
F2YP0BU02G7LK5 81 chr21 15557360 255 40M303H * 0 0 ATTTTATTCTCTTTGAAGCAATTGTGAATGGGAGTTTACT D>>>>IIIIIIHHG???IIIIIIIIIHHHFFEIH999HII
F2YP0BU02HXV58 145 chr21 15952091 255 226H40M6S * 0 0 GCAAATTGTGCTGCTTTAAACATGCGTGTGCAAGTATCTTtttcat AA===DDDDAAAAD???:::ABBBBBAAA:888ECF;F>>>?8??@
F2YP0BU02HREML 65 chr21 16386077 255 320H33M5S * 0 0 CCAAAGTTCTGGGATTACAGGCGTGAGCCATCGcgccc HH???HHIIIHFHIIIIIIICDDHHIIIIIIHHHHHHH
F2YP0BU02IOF1F 129 chr21 17567321 255 7S28M409H * 0 0 taaagagAAGAATTCTCAACCCAGAATTTCATATC 4100<<A>4113:<EFGGGFFFHHHHHHDFFFFED
F2YP0BU02IKX84 81 chr21 18491628 255 22M1D18M9S341H * 0 0 GTCTCTACCAAAAAATACAAAAATTAGCCGGGCGTGGTGGcatgtctgt ;;;=7@.55------?2?11112GGB=CCCCDIIIIIIIIIHHHHHHII
F2YP0BU02GW5VA 145 chr21 20255344 255 286H6S32M * 0 0 caagaaCAAACACATTCAAAAGCTAGTAGAAGGCAAGA IIIIIIIHHHIIIIIIICCCCIIIIIIIIIIIIIIIII
F2YP0BU02JIMJ4 65 chr21 22383051 255 208H19M * 0 0 CCCTTTATCATTTTTTATT 555544E?GE113344I22
F2YP0BU02IXZGF 145 chr21 23094798 255 291H13M1I18M * 0 0 GCAAGCTCCACTTCCCGGGTTCACGCCATTCT IIIIIIIIIIIGG;;;GGHIIIIIGGGIIIII
F2YP0BU02IODR5 129 chr21 30935325 255 37M154H * 0 0 GAAATAAAGGGTATTCAATTAGGAAAAGAGGAAGTCA 6...7/--..,30;9<<>@BFFFAAAAHIIIIIH@@@
F2YP0BU02IMZBL 145 chr21 31603486 255 342H28M1D1M * 0 0 ATACAAAAATTAGCCGGGCACAGTGGCAG BB1552222<<>9==8;;?AA=??A???A
F2YP0BU02JA9PR 145 chr21 31677159 255 229H23M * 0 0 CACACCTGTAACCCCAGCACTTT IIIIIIIIIIICCCCIIIIIHHH
F2YP0BU02HKC61 65 chr21 31678718 255 300H40M * 0 0 CACTGCACTCCAGCCTGGGTGACAAAGCAAGACTCTGTCT AA@BD:::==AAA@A?8888:<90004<>>?><<<<4442
F2YP0BU02HKC61 65 chr21 31678718 255 300H40M * 0 0 CACTGCACTCCAGCCTGGGTGACAAAGCAAGACTCTGTCT AA@BD:::==AAA@A?8888:<90004<>>?><<<<4442
F2YP0BU02HVA88 16 chr21 31703558 255 1M1D35M8S * 0 0 TGGGATTACAGGCGTGAGCTACCACACCCAGCCAGAgttcaaat >8888DFFHHGFHHHH@@?@?DDC96666HIIIFFFFFFFFFFF
F2YP0BU02JDCF1 129 chr21 31816600 255 38M103H * 0 0 AGGAGAATCGCTTGAACCCAGGAGGCAGAGGTTGCGGT IIIIIIIIIIIHHHIIHHHIIIIIIIIIIIIIIIIIII
F2YP0BU02GZ1GO 65 chr21 33360122 255 76H6S38M * 0 0 cctagaCTTCACACACACACACACACACACACACACACACACAC BBBBD?:688CFFFFFFFFFFFFFFFFFFFFFFFFFFDDBBB51
F2YP0BU02FX387 145 chr22 14786201 255 201H26M * 0 0 TGGATGAAGCTGGAAACCATCATTCT IIHHHHHHHHHHHHHFFFFFFFFFFF
F2YP0BU02IF2NE 65 chr22 16960842 255 209H40M10S * 0 0 TGGCATGCACCTGTAGTCTCAGCTACTTGGGAGGCTGAGGtgggaggatc BAAADDDDFDDDDDDBBA889<A?4444000@<>AA?9444;;8>77<7-
F2YP0BU02F4TVA 0 chr22 19200522 255 49M * 0 0 CCTGGGAGGCGGAGGTTGCAGTGAGCCGAGATCACGCCATTGCACTCCA FFF???FFFFFIIIIIIIIIIIIIIIIIIIIIIIHHIIFHFFFGDDB=5
F2YP0BU02HKC61 81 chr22 29516998 255 8S32M300H * 0 0 agacagagTCTTGCTTTGTCACCCAGGCTGGAGTGCAGTG 2444<<<<>?>><40009<:8888?A@AAA==:::DB@AA
F2YP0BU02FS4EM 65 chr22 30159364 255 339H29M * 0 0 CTCCTGCCTCAGCCTCCCGAGTAGTTGGG IIIIHHEIIIIHHHH??=DDHIIIIIDDD
F2YP0BU02G197P 129 chr22 32044496 255 40M10S258H * 0 0 TTGTTGGACATTTGGGTTGGTTCCAAGTCTTTGCTATTGTgaataatgcc IIIIIIIIIIHHHHHHIIIIIIIIIIIII;;;IIIIIIIIIIIIIIIIII
F2YP0BU02FIING 16 chr22 45959944 255 3M1I11M1I26M * 0 0 AGCTATGGTACTGGCTATGAAAGCAGACACATAGACCAATGG :::9:32267=:114244/...446==<<<?@?:9::::AAA
F2YP0BU02GUB9L 145 chr22 49198404 255 176H16M1I20M * 0 0 CACCACGCTCGGCTAATTTTTGTATTTTTAGTAGAGA IIIIIIIIIHAAC;<</////@4F5778;IIIIIIII
"""
left_suffix = "_L"
right_suffix = "_R"
# Read the mapping
mate_to_read_dict = {}
i = 0
for i, line in enumerate( file( tmp_mates_mapping_file_name, 'rb' ) ):
line = line.strip()
if not line.startswith( "#" ):
fields = line.split()
if len( fields ) != 4:
skip_line( "num_fields", i+1, line )
continue
mate_name, read_name, s_offset, e_offset = fields
if mate_name in mate_to_read_dict:
skip_line( 'two_mate_names', i+1, mate_name )
continue
mate_to_read_dict[ mate_name ] = ( read_name, int( s_offset ), int( e_offset ) )
# Read sequence data
read_to_nucs_dict = {}
seqs = 0
fasta_reader = FastaReader( file( input2, 'rb' ) )
while True:
seq = fasta_reader.next()
if not seq:
break
seqs += 1
seq_text_upper = seq.text.upper()
if seq.name in read_to_nucs_dict:
if seq_text_upper != read_to_nucs_dict[ seq.name ]:
skip_line( 'inconsistent_reads', seqs, seq.name )
continue
read_to_nucs_dict[ seq.name ] = seq_text_upper
# Read quality data
def quality_sequences( f ):
seq_name = None
seq_quals = None
line_number = 0
for line in f:
line_number += 1
line = line.strip()
if line.startswith( ">" ):
if seq_name != None:
yield ( seq_name, seq_quals, seq_line )
seq_name = sequence_name( line )
seq_line = line_number
seq_quals = []
elif seq_name is None:
skip_line( 'no_header', line_number, line )
continue
else:
seq_quals += [ int( q ) for q in line.split() ]
if seq_name is not None:
yield ( seq_name, seq_quals, seq_line )
def sequence_name( s ):
s = s[ 1: ].strip()
if not s:
return ""
else:
return s.split()[ 0 ]
read_to_quals_dict = {}
# TODO: should we use Dan's fastaNamedReader here?
for seq_name, quals, line_number in quality_sequences( file( input4 ) ):
quals = samify_phred_scores( quals )
if seq_name in read_to_quals_dict:
if quals != read_to_quals_dict[ seq_name ]:
skip_line( 'inconsistent_reads', line_number, seq_name )
continue
if len( quals ) != len( read_to_nucs_dict[ seq_name ] ):
skip_line( 'inconsistent_read_lengths', line_number, seq_name )
continue
read_to_quals_dict[ seq_name ] = quals
# process the SAM file
tmp_align_file_names = ' '.join( tmp_align_file_name_list )
combined_chrom_file_name = get_tmp_file_name( suffix='combined_chrom' )
command = 'cat %s | grep -v "^@" | sort -k 1 > %s' % ( tmp_align_file_names, combined_chrom_file_name )
run_command( command )
fout = file( output, 'w+b' )
has_non_header = False
i = 0
for i, line in enumerate( file( combined_chrom_file_name, 'rb' ) ):
line = line.strip()
if line.startswith( "@" ):
if has_non_header:
skip_line( 'sam_headers', i+1, line )
continue
fout.write( "%s\n" % line )
continue
has_non_header = True
fields = line.split()
num_fields = len( fields )
if num_fields < SAM_MIN_COLUMNS:
skip_line( 'sam_min_columns', i+1, line )
continue
# Set flags for mates
try:
flag = int( fields[ SAM_FLAG_COLUMN ] )
except ValueError:
skip_line( 'sam_flag', i+1, line )
continue
if not( flag & ( BAM_FPAIRED + BAM_FREAD1 + BAM_FREAD2 ) == 0 ):
skip_line( 'reads_paired', i+1, line )
continue
mate_name = fields[ SAM_QNAME_COLUMN ]
unmap_it = False
half = None
if mate_name.endswith( left_suffix ):
flag += BAM_FPAIRED + BAM_FREAD2
fields[ SAM_FLAG_COLUMN ] = "%d" % flag
unmap_it = True
half = "L"
elif mate_name.endswith( right_suffix ):
flag += BAM_FPAIRED + BAM_FREAD1
fields[ SAM_FLAG_COLUMN ] = "%d" % flag
unmap_it = True
half = "R"
on_plus_strand = ( flag & BAM_FREVERSE == 0 )
# Convert position from mate to read by adding clipping to cigar
if not unmap_it:
read_name = mate_name
else:
try:
read_name, s_offset, e_offset = mate_to_read_dict[ mate_name ]
except KeyError:
skip_line( 'missing_mate', i+1, mate_name )
continue
cigar = fields[ SAM_CIGAR_COLUMN ]
cigar_prefix = None
cigar_suffix = None
if half == "L":
if on_plus_strand:
if s_offset > 0:
cigar_prefix = ( s_offset, "S" )
if e_offset > 0:
cigar_suffix = ( e_offset, "H" )
else:
if e_offset > 0:
cigar_prefix = ( e_offset, "H" )
if s_offset > 0:
cigar_suffix = ( s_offset, "S" )
elif half == "R":
if on_plus_strand:
if s_offset > 0:
cigar_prefix = ( s_offset, "H" )
if e_offset > 0:
cigar_suffix = ( e_offset, "S" )
else:
if e_offset > 0:
cigar_prefix = ( e_offset, "S" )
if s_offset > 0:
cigar_suffix = ( s_offset, "H" )
else:
if on_plus_strand:
if s_offset > 0:
cigar_prefix = ( s_offset, "S" )
if e_offset > 0:
cigar_suffix = ( e_offset, "S" )
else:
if e_offset > 0:
cigar_prefix = ( e_offset, "S" )
if s_offset > 0:
cigar_suffix = ( s_offset, "S" )
if cigar_prefix != None:
count, op = cigar_prefix
cigar = prefix_cigar( "%d%s" % ( count, op ), cigar )
if op == "S":
refPos = int( fields[ SAM_POS_COLUMN ] ) - count
fields[ SAM_POS_COLUMN ] = "%d" % refPos
if cigar_suffix != None:
count, op = cigar_suffix
cigar = suffix_cigar( cigar,"%d%s" % ( count, op) )
fields[ SAM_QNAME_COLUMN ] = read_name
fields[ SAM_CIGAR_COLUMN ] = cigar
# Fetch sequence and quality values, and flip/clip them
if read_name not in read_to_nucs_dict:
skip_line( 'missing_seq', i+1, read_name )
continue
nucs = read_to_nucs_dict[ read_name ]
if not on_plus_strand:
nucs = reverse_complement( nucs )
quals = None
if read_to_quals_dict != None:
if read_name not in read_to_quals_dict:
skip_line( 'missing_quals', i+1, read_name )
continue
quals = read_to_quals_dict[ read_name ]
if not on_plus_strand:
quals = reverse_string( quals )
cigar = split_cigar( fields[ SAM_CIGAR_COLUMN ] )
nucs, quals = clip_for_cigar( cigar, nucs, quals )
fields[ SAM_SEQ_COLUMN ] = nucs
if quals != None:
fields[ SAM_QUAL_COLUMN ] = quals
# Output the line
fout.write( "%s\n" % "\t".join( fields ) )
fout.close()
def prefix_cigar( prefix, cigar ):
ix = 0
while cigar[ ix ].isdigit():
ix += 1
if cigar[ ix ] != prefix[ -1 ]:
return prefix + cigar
count = int( prefix[ :-1 ] ) + int( cigar[ :ix ] )
return "%d%s%s" % ( count, prefix[ -1 ], cigar[ ix+1: ] )
def suffix_cigar( cigar, suffix ):
if cigar[ -1 ] != suffix[ -1 ]:
return cigar + suffix
ix = len( cigar ) - 2
while cigar[ix].isdigit():
ix -= 1
ix += 1
count = int( cigar[ ix:-1 ] ) + int( suffix[ :-1 ] )
return "%s%d%s" % ( cigar[ :ix ], count, suffix[ -1 ] )
def split_cigar( text ):
fields = []
field = []
for ch in text:
if ch not in "MIDHS":
field += ch
continue
if field == []:
raise ValueError
fields += [ ( int( "".join( field ) ), ch ) ]
field = []
if field != []:
raise ValueError
return fields
def clip_for_cigar( cigar, nucs, quals ):
# Hard clip prefix
count, op = cigar[0]
if op == "H":
nucs = nucs[ count: ]
if quals != None:
quals = quals[ count: ]
count, op = cigar[ 1 ]
# Soft clip prefix
if op == "S":
nucs = nucs[ :count ].lower() + nucs[ count: ]
# Hard clip suffix
count,op = cigar[ -1 ]
if op == "H":
nucs = nucs[ :-count ]
if quals != None:
quals = quals[ :-count ]
count, op = cigar[ -2 ]
# Soft clip suffix
if op == "S":
nucs = nucs[ :-count ] + nucs[ -count: ].lower()
return nucs, quals
def samify_phred_scores( quals ):
"""
Convert a decimal list of phred base-quality scores to a sam quality string.
Note that if a quality is outside the dynamic range of sam's ability to
represent it, we clip the value to the max allowed. SAM quality scores
range from chr(33) to chr(126).
"""
if min( quals ) >= 0 and max( quals ) <= 126-33:
return "".join( [ chr( 33 + q ) for q in quals ] )
else:
return "".join( [ chr( max( 33, min( 126, 33+q ) ) ) for q in quals ] )
def reverse_complement( nucs ):
complementMap = maketrans( "ACGTacgt", "TGCAtgca" )
return nucs[ ::-1 ].translate( complementMap )
def reverse_string( s ):
return s[ ::-1 ]
def __main__():
# Parse command line
# input1: a reference genome ( 2bit or fasta )
# input2: a collection of 454 paired end reads ( a fasta file )
# input3: a linker sequence ( a very small fasta file )
# input4: a base quality score 454 file ( qual454 )
parser = optparse.OptionParser()
parser.add_option( '', '--ref_name', dest='ref_name', help='The reference name to change all output matches to' )
parser.add_option( '', '--ref_source', dest='ref_source', help='The reference is cached or from the history' )
parser.add_option( '', '--ref_sequences', dest='ref_sequences', help='Number of sequences in the reference dataset' )
parser.add_option( '', '--source_select', dest='source_select', help='Use pre-set or cached reference file' )
parser.add_option( '', '--input1', dest='input1', help='The name of the reference file if using history or reference base name if using cached' )
parser.add_option( '', '--input2', dest='input2', help='The 454 reads file to align' )
parser.add_option( '', '--input3', dest='input3', help='The sequencing linker file' )
parser.add_option( '', '--input4', dest='input4', help='The base quality score 454 file' )
parser.add_option( '', '--output', dest='output', help='The output file' )
parser.add_option( '', '--lastz_seqs_file_dir', dest='lastz_seqs_file_dir', help='Directory of local lastz_seqs.loc file' )
( options, args ) = parser.parse_args()
# output version # of tool
try:
tmp = tempfile.NamedTemporaryFile().name
tmp_stdout = open( tmp, 'wb' )
proc = subprocess.Popen( args='lastz -v', shell=True, stdout=tmp_stdout )
tmp_stdout.close()
returncode = proc.wait()
stdout = None
for line in open( tmp_stdout.name, 'rb' ):
if line.lower().find( 'version' ) >= 0:
stdout = line.strip()
break
if stdout:
sys.stdout.write( '%s\n' % stdout )
else:
raise Exception
except:
sys.stdout.write( 'Could not determine Lastz version\n' )
if options.ref_name:
ref_name = '[nickname=%s]' % options.ref_name
else:
ref_name = ''
if options.ref_source == 'history':
# Reference is a fasta dataset from the history
try:
# Ensure there is at least 1 sequence in the dataset ( this may not be necessary ).
error_msg = "The reference dataset is missing metadata, click the pencil icon in the history item and 'auto-detect' the metadata attributes."
ref_sequences = int( options.ref_sequences )
if ref_sequences < 1:
stop_err( error_msg )
except:
stop_err( error_msg )
else:
ref_sequences = 0
tmp_w12_name = get_tmp_file_name( suffix='vs_linker.W12' )
tmp_T1_name = get_tmp_file_name( suffix='vs_linker.T1' )
# Run lastz twice ( with different options ) on the linker sequence and paired end reads,
# looking for the linker ( each run finds some the other doesn't )
command = 'lastz %s %s W=12 --notrans --exact=18 --match=1,3 O=1 E=3 Y=10 L=18 --ambiguousn --coverage=85 --format=general-:name2,zstart2+,length2,size2 > %s' % \
( options.input3, options.input2, tmp_w12_name )
run_command( command )
command = 'lastz %s %s T=1 --match=1,2 O=1 E=2 X=15 K=10 Y=15 L=18 --ambiguousn --coverage=85 --format=general-:name2,zstart2+,length2,size2 > %s' % \
( options.input3, options.input2, tmp_T1_name )
run_command( command )
# Combine the alignment output from the two lastz runs
tmp_combined_linker_file_name = get_tmp_file_name( suffix='vs_linker' )
command = 'cat %s %s | sort -u > %s' % ( tmp_w12_name, tmp_T1_name, tmp_combined_linker_file_name )
run_command( command )
# Use the alignment info to split reads into left and right mates
tmp_mates_mapping_file_name, tmp_mates_file_name, tmp_mates_short_file_name, tmp_mates_long_file_name = split_paired_reads( options.input2, tmp_combined_linker_file_name )
# Align mates to the reference - tmp_align_file_names is a list of file names created by align_mates()
tmp_align_file_name_list = align_mates( options.input1, options.ref_source, ref_name, ref_sequences, tmp_mates_short_file_name, tmp_mates_long_file_name )
# Combine and convert mate coordinates back to read coordinates
paired_mate_unmapper( options.input2, options.input4, tmp_mates_mapping_file_name, tmp_align_file_name_list, options.output )
# Delete all temporary files
for file_name in tmp_file_names:
os.remove( file_name )
# Handle any invalid lines in the input data
if total_skipped_lines:
msgs = dict( bad_interval="Bad interval in line",
inconsistent_read_lengths="Inconsistent read/quality lengths for seq #",
inconsistent_reads="Inconsistent reads for seq #",
inconsistent_sizes="Inconsistent sizes for seq #",
missing_mate="Mapping file does not include mate on line",
missing_quals="Missing quality values for name on line",
missing_seq="Missing sequence for name on line",
multiple_seqs="Multiple names for seq #",
no_header="First quality sequence has no header",
num_fields="Must have 4 fields in line",
reads_paired="SAM flag indicates reads already paired on line",
sam_flag="Bad SAM flag on line",
sam_headers="SAM headers on line",
sam_min_columns="Need 11 columns on line",
two_mate_names="Mate name already seen, line",
wrong_seq_len="Size differs from length of seq #" )
print "Skipped %d invalid lines: "
msg = ""
for k, v in skipped_lines.items():
if v[0]:
# v[0] is the number of times the error occurred
# v[1] is the position of the line or sequence in the file
# v[2] is the name of the sequence or the text of the line
msg += "(%d)%s %d:%s. " % ( v[0], msgs[k], v[1], v[2] )
print msg
if __name__=="__main__": __main__()
Reference in New Issue View Git Blame Copy Permalink