Class: Transrate::Snap

Inherits:
Object
  • Object
show all
Defined in:
lib/transrate/snap.rb

Instance Attribute Summary collapse

Instance Method Summary collapse

Constructor Details

#initialize ⇒ Snap 



13
14
15
16
17
18
19
20
21
22
23
 
# File 'lib/transrate/snap.rb', line 13

def initialize
  which_snap = Cmd.new('which snap-aligner')
  which_snap.run
  if !which_snap.status.success?
    raise SnapError.new("could not find snap in the path")
  end
  @snap = which_snap.stdout.split("\n").first

  @index_built = false
  @index_name = ""
end

Instance Attribute Details

#bam ⇒ Object (readonly)

Returns the value of attribute bam.



11
12
13
 
# File 'lib/transrate/snap.rb', line 11

def bam
  @bam
end

#index_name ⇒ Object (readonly)

Returns the value of attribute index_name.



11
12
13
 
# File 'lib/transrate/snap.rb', line 11

def index_name
  @index_name
end

#read_count ⇒ Object (readonly)

Returns the value of attribute read_count.



11
12
13
 
# File 'lib/transrate/snap.rb', line 11

def read_count
  @read_count
end

Instance Method Details

#build_index(file, threads) ⇒ Object 



141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
 
# File 'lib/transrate/snap.rb', line 141

def build_index file, threads
  @index_name = File.basename(file, File.extname(file))
  unless Dir.exists?(@index_name)
    cmd = "#{@snap} index #{file} #{@index_name}"
    cmd << " -s 23"
    cmd << " -t#{threads}"
    cmd << " -bSpace" # contig name terminates with space char
    runner = Cmd.new cmd
    runner.run
    if !runner.status.success?
      err = runner.stderr
      msg = "Failed to build Snap index\n#{runner.stderr}"
      raise SnapError.new(msg)
    end
  end
  @index_built = true
end

#build_paired_cmd(l, r, threads) ⇒ Object 



25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
 
# File 'lib/transrate/snap.rb', line 25

def build_paired_cmd l, r, threads
  cmd = "#{@snap} paired #{@index_name}"
  l.split(",").zip(r.split(",")).each do |left, right|
    cmd << " #{left} #{right}"
  end
  cmd << " -o #{@bam}"
  cmd << " -s 0 1000" # min and max distance between paired-read starts
  cmd << " -H 300000" # max seed hits to consider in paired mode
  cmd << " -h 2000" # max seed hits to consider when reverting to single
  cmd << " -d 30" # max edit distance (function of read length?)
  cmd << " -t #{threads}"
  cmd << " -b" # bind threads to cores
  cmd << " -M"  # format cigar string
  cmd << " -D 5" # extra edit distance to search. needed for -om
  cmd << " -om 5" # Output multiple alignments. extra edit distance
  cmd << " -omax 10" # max alignments per pair/read
  cmd
end

#load_readcount(reads) ⇒ Object 



120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
 
# File 'lib/transrate/snap.rb', line 120

def load_readcount reads
  @read_count = 0
  if File.exist?("#{@read_count_file}")
    @read_count = File.open("#{@read_count_file}").readlines.join.to_i
  else
    reads.split(",").each do |l|
      cmd = "wc -l #{l}"
      count = Cmd.new(cmd)
      count.run
      if count.status.success?
        @read_count += count.stdout.strip.split(/\s+/).first.to_i/4
        File.open("#{@read_count_file}", "wb") do |out|
          out.write("#{@read_count}\n")
        end
      else
        logger.warn "couldn't get number of reads from #{l}"
      end
    end
  end
end

#map_reads(file, left, right, outputname: nil, threads: 8) ⇒ Object

Raises:



44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
# File 'lib/transrate/snap.rb', line 44

def map_reads(file, left, right, outputname: nil, threads: 8)
  raise SnapError.new("Index not built") if !@index_built

  lbase = File.basename(left.split(",").first)
  rbase = File.basename(right.split(",").first)
  index = File.basename(@index_name)
  @bam = File.expand_path("#{lbase}.#{rbase}.#{index}.bam")
  @read_count_file = "#{lbase}-#{rbase}-read_count.txt"

  @fixer = Fixer.new # from the fix-trinity-output gem
  unless File.exists? @bam
    snapcmd = build_paired_cmd(left, right, threads)
    runner = Cmd.new snapcmd
    runner.run
    save_readcount runner.stdout
    save_logs(runner.stdout, runner.stderr)
    unless runner.status.success?
      if runner.stderr=~/Unmatched\sread\sIDs/
        logger.warn runner.stderr
        logger.warn "Unmatched read IDs. Fixing input files..."
        remap_reads(left, right, threads)
      else
        raise SnapError.new("Snap failed\n#{runner.stderr}")
      end
    end
  else
    load_readcount left
  end
  @bam
end

#remap_reads(left, right, threads) ⇒ Object 



83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
 
# File 'lib/transrate/snap.rb', line 83

def remap_reads(left, right, threads)
  fixedleft = []
  fixedright = []
  i = 0
  left.split(",").zip(right.split(",")).each do |l, r|
    prefix = "reads-#{i}"
    @fixer.run(l, r, "#{prefix}")
    fixedleft << "#{prefix}-fixed.1.fastq"
    fixedright << "#{prefix}-fixed.2.fastq"
    i+=1
  end
  left = fixedleft.join(",")
  right = fixedright.join(",")
  File.delete(@bam)
  logger.info "Fixed input files"
  snapcmd = build_paired_cmd(left, right, threads)
  runner = Cmd.new snapcmd
  runner.run
  save_readcount runner.stdout
  save_logs(runner.stdout, runner.stderr)
  unless runner.status.success?
    raise SnapError.new("Snap failed\n#{runner.stderr}")
  end
end

#save_logs(stdout, stderr) ⇒ Object 



75
76
77
78
79
80
81
 
# File 'lib/transrate/snap.rb', line 75

def save_logs(stdout, stderr)
  FileUtils.mkdir_p 'logs'
  File.open('logs/snap.log', 'a') do |f|
    f.write stdout
    f.write stderr
  end
end

#save_readcount(stdout) ⇒ Object 



108
109
110
111
112
113
114
115
116
117
118
 
# File 'lib/transrate/snap.rb', line 108

def save_readcount stdout
  stdout.split("\n").each do |line|
    cols = line.split(/\s+/)
    if cols.size > 5 and cols[0]=~/[0-9\,]+/
      @read_count = cols[0].gsub(",", "").to_i / 2
      File.open("#{@read_count_file}", "wb") do |out|
        out.write("#{@read_count}\n")
      end
    end
  end
end