Skip to content

Commit

Permalink
No Idea what I changed here.
Browse files Browse the repository at this point in the history
  • Loading branch information
stela2502 committed Oct 24, 2024
1 parent a7cd883 commit 09e7ff9
Show file tree
Hide file tree
Showing 5 changed files with 16 additions and 82 deletions.
2 changes: 1 addition & 1 deletion src/analysis/gene_mapper.rs
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,7 @@ use std::thread;
use rayon::prelude::*;
use rayon::slice::ParallelSlice;

use crate::ofiles::{Ofiles, Fspot};
use crate::ofiles::{Ofiles};
//use std::fs::write;

//static EMPTY_VEC: Vec<String> = Vec::new();
Expand Down
2 changes: 1 addition & 1 deletion src/analysis/genomic_mapper.rs
Original file line number Diff line number Diff line change
Expand Up @@ -44,7 +44,7 @@ use std::thread;
use rayon::prelude::*;
use rayon::slice::ParallelSlice;

use crate::ofiles::{Ofiles, Fspot};
use crate::ofiles::{Ofiles};
//use std::fs::write;
use glob::glob;

Expand Down
22 changes: 10 additions & 12 deletions src/gene.rs
Original file line number Diff line number Diff line change
Expand Up @@ -76,14 +76,12 @@ impl Gene{
Err(e) => panic!("I could not parse the end of the transcript as usize: {e:?}"),
};

let sense_strand = sense_strand_s == "+";

Self{
chrom: chrom.to_string(),
start,
end,
exons,
sense_strand,
sense_strand: match sense_strand_s{ "+" => true, _ => false },
name: name.to_string(),
transcript: transcript.to_string(),
ids,
Expand Down Expand Up @@ -165,10 +163,10 @@ impl Gene{

/// generate mRNA and nacent RNA (if the last exon is small enough)
pub fn generate_rna_and_nascent_strings(&self, seq: &[u8], covered_area: usize) -> (Option<Vec<u8>>, Option<Vec<u8>>) {
let mrna = self.to_mrna(seq.to_owned(), covered_area);
let nascent = self.to_nascent( seq, covered_area );
let mrna = self.to_mrna( seq, covered_area );

let nascent = self.to_nascent(seq.to_owned(), covered_area);
(mrna, nascent)
( mrna, nascent)
}

/// returns the start position on the mRNA in genomic coordinates and the global end position
Expand Down Expand Up @@ -196,7 +194,7 @@ impl Gene{
/// get the mRNA sequence of the transcript in sense orientation.
/// Fails if any other base than AGCT is in the sequence
/// This returns the revers complement if on the opposite starnd
pub fn to_mrna(&self, seq: Vec<u8>, covered_area:usize) -> Option<Vec<u8>> {
pub fn to_mrna(&self, seq: &[u8], covered_area:usize) -> Option<Vec<u8>> {

let mut mrna = Vec::<u8>::with_capacity(seq.len() ); // Allocate more space for potential additions

Expand All @@ -216,12 +214,12 @@ impl Gene{
mrna.extend_from_slice(&seq[reg[0] - 1..(reg[1])]);
}
//println!("the final mRNA {:?}", std::str::from_utf8(&mrna) );
self.cut_to_size(mrna, covered_area )
self.cut_to_size( &mrna, covered_area )

}

/// cut the RNA to the right size
fn cut_to_size( &self, seq:Vec<u8>, covered_area:usize) -> Option<Vec<u8>> {
fn cut_to_size( &self, seq: &[u8], covered_area:usize) -> Option<Vec<u8>> {

let start = seq.len().saturating_sub(covered_area);
if ! self.sense_strand{
Expand All @@ -244,7 +242,7 @@ impl Gene{
/// get the nascent RNA for this transcript (including introns).
/// Fails if any other base than AGCT is in the sequence
/// This returns the revers complement if on the opposite starnd
fn to_nascent(&self, seq:Vec<u8> , covered_area:usize) -> Option<Vec<u8>> {
fn to_nascent(&self, seq: &[u8] , covered_area:usize) -> Option<Vec<u8>> {
if self.exons.len() == 0{
eprintln!("I have no exons?! - something is wrong here! {self}");
return None
Expand All @@ -266,15 +264,15 @@ impl Gene{
let (_glob_start, _glob_end) = self.to_mrna_positions(covered_area);
//make the nascent WAY longer??
//self.cut_to_size(nascent, glob_end - glob_start )
self.cut_to_size(nascent, covered_area*2 )
self.cut_to_size(&nascent, covered_area*2 )
}
else {
None
}
}

/// the reverse complement of a Vec<u8>
fn rev_compl( seq:Vec<u8> ) -> Vec<u8>{
fn rev_compl( seq:&[u8] ) -> Vec<u8>{
seq.iter()
.rev()
.filter_map(|&c| COMPLEMENT[c as usize])
Expand Down
70 changes: 3 additions & 67 deletions src/genes_mapper/genes_mapper.rs
Original file line number Diff line number Diff line change
Expand Up @@ -416,67 +416,6 @@ impl GenesMapper{
//println!("\nMakes sense? {:?}", data);
data
}

/*fn query_database( &self, seq: &[u8] ) -> Vec<((usize, i32), usize)> {
//if self.small_entries {
return self.query_database_small( seq );
//}
let mut res = HashMap::<(usize, i32), usize>::new();// Vec<i32>>::new();
let mut read_data = GeneData::new( seq,"read_t", "read", "read2", 0 );
let mut i = 0;
let mut failed = 0;
while let Some((key, start)) = read_data.next(){
//println!("after {i} (matching) iterations ({}) the res looks like that {res:?}", Self::key_to_string( &key) );
if self.reject_key(&key){
continue;
}
if ! self.mapper[key as usize].is_empty() {
i+=1;
failed = 0;
#[cfg(debug_assertions)]
if self.debug {
let values: Vec<_> = self.mapper[key as usize].data().collect();
let m = 9.min( values.len() );
println!("I am testing this key: {} and got {:?} ...", Self::as_dna_string(&key) , &values[0..m] );
}
// this will directly modify the res HashMap as it can also find more than one!
self.mapper[key as usize].get( &mut res, start as i32 );
}
else {
failed +=1;
if failed == 3 {
match read_data.iterator_skip_to_next_frame() {
Ok(_) => failed = 0,
Err(_) => break,
}
}
#[cfg(debug_assertions)]
if self.debug{
println!("I am testing this key: {} and got Nothing", Self::as_dna_string(&key) );
}
}
if i == 30 {
break;
}
}
let mut res_vec: Vec<_> = res.into_iter().collect();
res_vec.retain(|(_, counts)| *counts > 4);
// longest first
res_vec.sort_by(|( (a_key, _a_rel_start), a_counts), ((b_key, _b_rel_start), b_counts)| {
match b_counts.cmp(&a_counts) {
std::cmp::Ordering::Equal => a_key.cmp(b_key),
other => other,
}
});
res_vec
}
fn query_database_small( &self, seq: &[u8] ) -> Vec<((usize, i32), usize)> {
*/

fn query_database( &self, seq: &[u8] ) -> Vec<((usize, i32), usize)> {
let mut res = HashMap::<(usize, i32), usize>::new();// Vec<i32>>::new();
Expand Down Expand Up @@ -618,16 +557,13 @@ impl GenesMapper{
break
}
}
#[cfg(debug_assertions)]
if self.debug{
println!("I got this nw: {nw} and the matches up to now: {}",helper);
}
};

} // end populating helper

#[cfg(debug_assertions)]
if self.debug{
//let res_lines = res_vec.iter().map(|obj| format!("{:?}", obj)).collect::<Vec<String>>().join("\n");
{
//let res_lines = res_vec.iter().map(|obj| format!("{:?}", obj)).collect::<Vec<String>>().join("\n");
println!("mapping the read {read_data}\nwe obtained an initial mapping result \n{helper}\nand in the end found \n{:?}",
helper.get_best( seq.len()) );
println!("#################################################################################");
Expand Down
2 changes: 1 addition & 1 deletion src/genes_mapper/needleman_wunsch_affine.rs
Original file line number Diff line number Diff line change
Expand Up @@ -383,7 +383,7 @@ impl <'a> NeedlemanWunschAffine {
let mut drop_replaces = 0;

// Check nucleotide matches and adjust Cigar vector
while let (Some(nuc1), Some(nuc2)) = (read.get_nucleotide_2bit(read_id.saturating_sub(1)), database.get_nucleotide_2bit(database_id.saturating_sub(1) )) {
while let (Some(_nuc1), Some(_nuc2)) = (read.get_nucleotide_2bit(read_id.saturating_sub(1)), database.get_nucleotide_2bit(database_id.saturating_sub(1) )) {

matching += 1;

Expand Down

0 comments on commit 09e7ff9

Please sign in to comment.