Class: Mapping

Inherits:
Object
  • Object
show all
Defined in:
lib/full_lengther_next/mapping.rb

Instance Method Summary collapse

Constructor Details

#initialize(user_options = {}) ⇒ Mapping

Returns a new instance of Mapping.



57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
 
# File 'lib/full_lengther_next/mapping.rb', line 57

def initialize(user_options = {})

	options = {
		threads: 1,
		total_reads: [],
	}.merge!(user_options)
	@ref_fasta_path = options[:ref_fasta_path]
	@temp_folder = options[:temp_folder]
	@threads = options[:threads]

	@map_files = []
	@paired = []
	@idxstats = []
	@mpileups = []
	@coverage_results = {}
	@total_reads = options[:total_reads]
end

Instance Method Details

#do_map(user_options = {}) ⇒ Object 



96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
 
# File 'lib/full_lengther_next/mapping.rb', line 96

def do_map(user_options = {})
	options = {
		files: [],
		command: 'bowtie2 -p /THREADS/ -x /REFERENCE/',
		paired_pipe: '| samtools view -bS -F 4 | samtools sort -o /OUTPUT/',
		single_pipe: '| samtools view -bS -F 4 | samtools sort -o /OUTPUT/',
		flag_single: '-U',
		flags_paired: ['-1', '-2'],
		additional_paired_flags: '',
		flag_output: '-S',
		output: File.join(@temp_folder, 'map_data'),
		log: File.join(@temp_folder, 'mapping_log'),
		force: false
	}
	options.merge!(user_options)
	options[:files].each_with_index do |read_files, map_process_id|
		cmd = options[:command].gsub('/THREADS/', @threads.to_s)
		cmd.gsub!('/REFERENCE/', @ref)
		if read_files.length == 1
			cmd = cmd + " #{options[:flag_single]} #{read_files.first}"
			@paired << false
		elsif read_files.length == 2
			@paired << true
			cmd = cmd + " #{options[:additional_paired_flags]} #{options[:flags_paired].first} #{read_files.first} #{options[:flags_paired].last} #{read_files.last}" 
		else
			raise('Incorrect number of read files. Must be 1 (single) or 2 (paired).')
		end
		map_file = nil
		if options[:paired_pipe].nil? || options[:single_pipe].nil?
			map_file = options[:output] + "_#{map_process_id}" + '.sam'
			cmd = cmd + " #{options[:flag_output]} #{map_file} &> #{options[:log]}_#{map_process_id}"
		else
			if @paired[map_process_id]
				pipe = options[:paired_pipe]
			else
				pipe = options[:single_pipe]
			end
			map_file = options[:output] + "_#{map_process_id}" + '.bam'
			cmd = cmd + " 2> #{options[:log]}_#{map_process_id} " + pipe.gsub('/OUTPUT/', map_file)
		end
		@map_files << map_file
		system(cmd) if options[:force] || !File.exists?(map_file)
		@total_reads << File.open("#{options[:log]}_#{map_process_id}").readlines.select{|line| /\d+ reads; of these:/ =~ line}.first.split(' ').first.to_i if File.exists?("#{options[:log]}_#{map_process_id}") && @total_reads[map_process_id].nil?
		raise('ERROR: The mapping process has failed, please check the map folder into the temp folder') if @total_reads[map_process_id].nil? || @total_reads[map_process_id] == 0
	end
end

#do_ref(user_options = {}) ⇒ Object 



75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
 
# File 'lib/full_lengther_next/mapping.rb', line 75

def do_ref(user_options = {})
	options = {
		name: 'ref',
		command: 'bowtie2-build -f --threads /THREADS/ /REF_FASTA/ /OUTPUT/',
		log: 'reference_log',
		force: false
	}
	options.merge!(user_options)
	@ref = File.join(@temp_folder, options[:name])
	cmd = options[:command].gsub('/THREADS/', @threads.to_s)
	cmd.gsub!('/REF_FASTA/', @ref_fasta_path)
	cmd.gsub!('/OUTPUT/', @ref)
	cmd = cmd + " &> #{File.join(@temp_folder, options[:log])}"
	system(cmd) if options[:force] || Dir.glob(@ref+'*.bt2').length == 0 
end

#do_samtools_ref ⇒ Object 



91
92
93
94
 
# File 'lib/full_lengther_next/mapping.rb', line 91

def do_samtools_ref
	cmd = "samtools faidx #{@ref_fasta_path}"
	system(cmd) if !File.exists?(@ref_fasta_path + '.fai')
end

#idxstats ⇒ Object 



159
160
161
162
163
164
165
166
167
 
# File 'lib/full_lengther_next/mapping.rb', line 159

def idxstats
	@map_files.each_with_index do |map_file, map_process_id|
		prefix = File.basename(map_file).gsub(/\.bam|\.sam|\.cram/, '')
		file_path = File.join(@temp_folder, "#{prefix}_idxstats_#{map_process_id}.gz")
		cmd = "samtools idxstats #{map_file} | gzip - -f > #{file_path}"
		system(cmd) if !File.exists?(file_path)
		parse_idxstats(file_path)
	end
end

#index(user_options = {}) ⇒ Object 



143
144
145
146
147
 
# File 'lib/full_lengther_next/mapping.rb', line 143

def index(user_options = {})
	@map_files.each do |map_file|
		system("samtools index #{map_file}") if (map_file.include?('.bam') && !File.exists?(map_file+'.bai')) || user_options[:force]
	end
end

#mpileup(user_options = {}) ⇒ Object 



169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
 
# File 'lib/full_lengther_next/mapping.rb', line 169

def mpileup(user_options = {})
	parse_options = {
		add_coverages: false, 
		normalize_coverages: false,
		cols: [1,2,4] # 1 based for cut
	}
	parse_options.merge!(user_options.delete(:parse_options)) if !user_options[:parse_options].nil?
	opts = []
	do_samtools_ref
	user_options.each do |flag, value|
		opts << [flag, value.to_s]
	end

	contig_list_file = File.join(@temp_folder, File.basename(@ref_fasta_path)+'.lst')
	system("grep '>' #{@ref_fasta_path} | sed 's/>//g' > #{contig_list_file}") if !File.exists?(contig_list_file)
	idxstats if @idxstats.empty?
	cut = nil
	cut = " |cut -f #{parse_options[:cols].join(',')}" if !parse_options[:cols].nil? && !parse_options[:cols].empty?
	mpileup_files = []
	@map_files.each_with_index do |map_file, map_process_id|
		prefix = File.basename(map_file).gsub(/\.bam|\.sam|\.cram/, '')
		file_path = File.join(@temp_folder, "#{prefix}_mpileup_#{map_process_id}.gz")
		mpileup_files << file_path
		cmd = "samtools mpileup -f #{@ref_fasta_path} #{opts.join(' ')} #{map_file}#{cut} | gzip - -f > #{file_path}" 
		system(cmd) if !File.exists?(file_path)
	end
	coverage_results = {}

	parse_mpileup(mpileup_files, contig_list_file) do |contig_name, contig_length, coverages|
		mapped_reads = @idxstats.map{|info| info[contig_name][:mapped]}.inject { |sum, n| sum + n }
		get_coverage_parameters(contig_name, contig_length, mapped_reads, coverages, parse_options, coverage_results)
	end
	return coverage_results
end

#parse_idxstats(file_path) ⇒ Object 



221
222
223
224
225
226
227
228
229
230
231
 
# File 'lib/full_lengther_next/mapping.rb', line 221

def parse_idxstats(file_path)
	stats = {}
	stats_file = ScbiZcatFile.new(file_path)
	while !stats_file.eof
        fields = stats_file.readline.chomp.split("\t")
        stats[fields[0]] = {length: fields[1].to_i, mapped: fields[2].to_i, unmmapped: fields[3].to_i}
    end
    stats_file.close
    stats.delete('*')
    @idxstats << stats
end

#parse_mpileup(file_paths, contig_list_file) ⇒ Object 



204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
 
# File 'lib/full_lengther_next/mapping.rb', line 204

def parse_mpileup(file_paths, contig_list_file)
	last_contig = nil
	mpileup_files = file_paths.map{|file_path| Mpileup.new(file_path)}
	File.open(contig_list_file).each do |contig_name|
		contig_name.chomp!
		contig_length = @idxstats.first[contig_name][:length]
		all_coverages = []
		mpileup_files.each do |mpileup_file|
			coverages = mpileup_file.read_contig(contig_name, contig_length)
			all_coverages << coverages if !coverages.nil? && !coverages.empty?  
		end
		yield(contig_name, contig_length, all_coverages)
	end
	mpileup_files.map{|mf| mf.close}
end

#report ⇒ Object 



149
150
151
152
153
154
155
156
157
 
# File 'lib/full_lengther_next/mapping.rb', line 149

def report
	reports = []
	@map_files.each do |map_file|
		cmd = "samtools flagstat #{map_file}"
		report = %x[#{cmd}].split("\n")
		reports << report
	end
	return reports
end