Module: BioPangenome

Defined in:
lib/bio-pangenome/pangenome.rb

Defined Under Namespace

Classes: GeneFlankingRegion, Transcript

Class Method Summary collapse

Class Method Details

.align_gene_groups(seqs: {}, tmp_folder: "/Volumes/PanGenome/GeneRegions/201910_v2_v3/tmp", output: "../pairwise_blast_oct_2019/varieties_6A_identites", distance: 0) ⇒ Object 



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/bio-pangenome/pangenome.rb', line 103

def self.align_gene_groups( seqs:{}, tmp_folder:"/Volumes/PanGenome/GeneRegions/201910_v2_v3/tmp", output:"../pairwise_blast_oct_2019/varieties_6A_identites", distance: 0 )
  out_tmp="#{tmp_folder}/out.blast"
  FileUtils.mkdir_p(tmp_folder)
  out = File.open("#{output}_#{distance}bp.tab", "w")
  out.puts [ "transcript" , "query", "subject" ,  "var_query", "var_subject", "aln_type", "length" , "pident" , "Ns_query", "Ns_subject", "Ns_total", "Flanking"   ].join("\t")
  seqs.each_pair do |transcript, transcript_seqs|
    vars = transcript_seqs.keys
    vars_done = []
    ns = {}
    vars.each do |v1|
      tmp =  tmp_folder  + "/" + v1 + ".fa"
      s = transcript_seqs[v1]
      seq = ">#{s.id}\n#{s.sequence}"
      File.open(tmp, 'w') {|f| f.write(seq) }
      ns[v1] = s.sequence.count('Nn')
    end
    vars.each do |v1|
      tmp1 =  tmp_folder  + "/" + v1 + ".fa"
      s1 = transcript_seqs[v1]
      next unless s1.sequence.length > 0
      vars.each do |v2|
        next if v1 == v2
        next if vars_done.include? v2
        s2 = transcript_seqs[v2]
        next unless s2.sequence.length > 0
        tmp2 =  tmp_folder  + "/" + v2 + ".fa"
        to_print = [transcript, s1.id , s2.id , v1,v2,"#{v1}->#{v2}"]
        to_print << blast_pair_fast(tmp1, tmp2, out_tmp) 
        to_print << ns[v1] 
        to_print << ns[v2]
        to_print << ns[v1] + ns[v2]
        to_print << distance
        out.puts to_print.join("\t")
      end
      vars_done << v1
    end
  end
  out.close
end

.blast_pair_fast(path_a, path_b, out_path, program: "blastn") ⇒ Object 



58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
 
# File 'lib/bio-pangenome/pangenome.rb', line 58

def self.blast_pair_fast(path_a, path_b, out_path, program: "blastn")
  cmd = "#{program} -query #{path_a} -subject #{path_b} -task #{program} -out #{out_path} -outfmt '5' "
  system cmd
  n = Bio::BlastXMLParser::XmlIterator.new(out_path).to_enum
  max_length = 0
  max_pident = 0.0
  n.each do | iter |
    iter.each do | hit |
      hit.each do | hsp |
        if hsp.align_len > max_length
          max_length = hsp.align_len
          max_pident = 100 * hsp.identity.to_f / hsp.align_len.to_f
        end
      end
    end
  end
  [max_length, max_pident]
end

.load_cds_sequences(varieties: [], genes: {}, prefix: "../flanking/filtered/", suffix: ".cds.fa.gz", set_id: "cds") ⇒ Object 



143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
 
# File 'lib/bio-pangenome/pangenome.rb', line 143

def self.load_cds_sequences( varieties:[], genes:{}, prefix: "../flanking/filtered/",  suffix: ".cds.fa.gz", set_id: "cds" )
  ret = Hash.new { |h, k| h[k] = Hash.new }
  varieties.each do |variety|
    path = "#{prefix}/#{variety}#{suffix}"
    infile = open(path)
    io = Zlib::GzipReader.new(infile) 
    Bio::FlatFile.open(Bio::FastaFormat, io) do |fasta_file|
      fasta_file.each do |entry|
        arr = entry.definition.split(".")
        next unless genes[arr[0]]
        row = genes[arr[0]]
        seq_name = GeneFlankingRegion.new(entry.definition,
          row["gene"], "",
          "", entry.definition, set_id, nil, variety )
        seq = entry.seq
        seq.gsub!(/N*$/, '')
        seq.gsub!(/^N*/, '')
        seq_name.sequence = seq
        base_gene = seq_name.gene
        ret[base_gene][variety] = seq_name unless ret[base_gene][variety]
      end
    end
    io.close
  end
  ret
end

.load_genes(filename, window: 0, no_windows: 0) ⇒ Object 



182
183
184
185
186
187
188
189
190
191
192
193
 
# File 'lib/bio-pangenome/pangenome.rb', line 182

def self.load_genes(filename, window: 0, no_windows: 0)
  genes = File.readlines(filename).map do |t|
    t.chomp!.split(".")[0]
  end
  if no_windows > 0
    puts "'loading window #{window} of #{no_windows}'"
    window_size = genes.size/no_windows
    start = window * window_size
    genes = genes[start, window_size]
  end
  genes
end

.load_lines(filename) ⇒ Object 



195
196
197
198
199
 
# File 'lib/bio-pangenome/pangenome.rb', line 195

def self.load_lines(filename)
  File.readlines(filename).map do |t|
    t.chomp!.rstrip
  end
end

.load_mapping_hash(varieties: [], transcripts: [], genes: [], distance: 1000, prefix: "../flanking/releasePGSBV1/", suffix: ".RefSeqv1.1") ⇒ Object 



38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
 
# File 'lib/bio-pangenome/pangenome.rb', line 38

def self.load_mapping_hash(varieties:[],  transcripts:[], genes:[], distance: 1000, prefix: "../flanking/releasePGSBV1/",  suffix: ".RefSeqv1.1")
  ret = Hash.new { |h, k| h[k] = Hash.new }
  varieties.each do |v|
    path = "#{prefix}#{distance}bp/#{v}_#{distance}bp_#{suffix}.reg.map"
    $stderr.puts path
    File.foreach(path) do |line|
      line.chomp!
      arr = line.split("\t")
      begin
        parsed = parseSequenceName(arr[0], arr[1])
      rescue Exception => e
        throw "Unable to parse #{line} (#{v}) [#{e.to_s}]" 
      end
      next unless transcripts.include? parsed.transcript or genes.include? parsed.gene
      ret[v][parsed.region] = parsed
    end
  end
  ret
end

.load_projected_genes(transcript_mapping, genes: []) ⇒ Object 



170
171
172
173
174
175
176
177
178
179
180
 
# File 'lib/bio-pangenome/pangenome.rb', line 170

def self.load_projected_genes(transcript_mapping, genes:[])
  projected_genes = {}
  Zlib::GzipReader.open(transcript_mapping) do |gzip|
    csv = CSV.new(gzip, headers: true)
    csv.each do |row|
      next unless genes.include? row["gene"]
      projected_genes[row["projected_gene"]] = row
    end
  end
  projected_genes
end

.load_sequences_from_hash(coordinates: {}, prefix: "../flanking/filtered/", suffix: "RefSeqv1.1", distance: 1000, projected_genes: {}) ⇒ Object 



78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
 
# File 'lib/bio-pangenome/pangenome.rb', line 78

def self.load_sequences_from_hash(coordinates:{},  prefix: "../flanking/filtered/",  suffix: "RefSeqv1.1", distance: 1000, projected_genes: {})
  ret = Hash.new { |h, k| h[k] = Hash.new }
  coordinates.each_pair do |variety, coords|

    path = "#{prefix}/#{distance}bp/#{variety}_#{distance}bp_#{suffix}.fa.gz"
    puts "Loading: #{path}"
    infile = open(path)
    io = Zlib::GzipReader.new(infile) 
    Bio::FlatFile.open(Bio::FastaFormat, io) do |fasta_file|
      fasta_file.each do |entry|
        next unless coords[entry.definition]
        seq_name = coords[entry.definition]
        seq = entry.seq
        seq.gsub!(/N*$/, '')
        seq.gsub!(/^N*/, '')
        seq_name.sequence = seq
        base_gene = projected_genes[seq_name.gene]["gene"]
        ret[base_gene][variety] = seq_name unless ret[base_gene][variety]
      end
    end
    io.close
  end
  ret
end

.parseEITranscript(name) ⇒ Object 



18
19
20
21
22
23
 
# File 'lib/bio-pangenome/pangenome.rb', line 18

def self.parseEITranscript name
  arr=name.split(".")
  match = /Traes(?<chr>[[:upper:]]{3}_scaffold_[[:digit:]]*)_(?<ver>[[:digit:]]{2})G(?<count>[[:digit:]]+)(?<conf>[[:upper:]]*)/.match arr[0]
  raise "Unable to parse: #{name}" unless match
  Transcript.new(name, arr[0],match[:chr].downcase,match[:ver],match[:count],arr[1],match[:conf], match[:count].to_i, arr[1])
end

.parsePGSBTranscript(name) ⇒ Object 



25
26
27
28
29
30
31
 
# File 'lib/bio-pangenome/pangenome.rb', line 25

def self.parsePGSBTranscript name
  arr=name.split(".")
  match = /Traes(?<variety>[[:upper:]]{3})(?<chr>[[:alnum:]]{1,2})(?<ver>[[:digit:]]{2})G(?<count>[[:digit:]]+)(?<conf>[[:upper:]]*)/.match arr[0]

  raise "Unable to parse: #{name}" unless match
  Transcript.new(name, arr[0],match[:chr],match[:ver],match[:count],match[:variety], match[:conf],match[:count].to_i, arr[1])
end

.parseSequenceName(region, name) ⇒ Object 



32
33
34
35
36
 
# File 'lib/bio-pangenome/pangenome.rb', line 32

def self.parseSequenceName region, name
  match = /(?<transcript>[[:alnum:]].+)_(?<ann>.+)_(?<flank_length>[[:digit:]]+bp)/.match name
  arr2=match[:transcript].split "."
  GeneFlankingRegion.new(match[:transcript],arr2[0],match[:ann], region, name, match[:flank_length] , nil, nil)
end

.parseTranscript(name) ⇒ Object 



12
13
14
15
16
17
 
# File 'lib/bio-pangenome/pangenome.rb', line 12

def self.parseTranscript name
  arr=name.split(".")
  match = /TraesCS(?<chr>[[:alnum:]]{1,2})(?<ver>[[:digit:]]{2})G(?<count>[[:digit:]]+)(?<conf>[[:upper:]]*)/.match arr[0]
  raise "Unable to parse: #{name}" unless match
  Transcript.new(name, arr[0],match[:chr],match[:ver],match[:count],arr[1],match[:conf], match[:count].to_i, arr[1])
end