Skip to content

Commit

Permalink
update manuscript demo code to use rust, keep all fields in the lifto…
Browse files Browse the repository at this point in the history
…ver output
  • Loading branch information
pre-mRNA committed Apr 11, 2024
1 parent ed3826d commit 1c4e079
Show file tree
Hide file tree
Showing 3 changed files with 9 additions and 9 deletions.
4 changes: 2 additions & 2 deletions manuscript_demo/process_data.sh
Original file line number Diff line number Diff line change
Expand Up @@ -17,10 +17,10 @@ export wd="/g/data/lf10/as7425/R2DTool_demo/"; mkdir -p ${wd}
bash ~/R2Dtool/scripts/cheui_to_bed.sh ${methylation_calls} "${wd}/methylation_calls.bed"

# annotate the methylation calls against gencode v38
r2d annotate -i "${wd}/methylation_calls.bed" -g ${annotation} -H > "${wd}/methylation_calls_annotated.bed"
time r2d annotate -i "${wd}/methylation_calls.bed" -g ${annotation} -H > "${wd}/methylation_calls_annotated.bed"

# liftover the annotated called to genomic coordinates
r2d liftover -i "${wd}/methylation_calls.bed" -g ${annotation} -H > "${wd}/methylation_calls_annotated.bed"
time r2d liftover -i "${wd}/methylation_calls_annotated.bed" -g ${annotation} -H > "${wd}/methylation_calls_annotated_lifted.bed"

# annotate methylation calls using Gencode v38
time Rscript ~/R2Dtool/scripts/R2_annotate.R "${wd}/methylationCalls.bed" "${annotation}" "${wd}/methylationCalls_annotated.bed"
Expand Down
4 changes: 2 additions & 2 deletions src/liftover.rs
Original file line number Diff line number Diff line change
Expand Up @@ -58,8 +58,8 @@ pub fn convert_transcriptomic_to_genomic_coordinates(
// Get the genomic strand
let genomic_strand = &transcript.strand;

// Join the additional columns (fields after the second one) using a tab character
let additional_columns = site_fields[2..].join("\t");
// Keep all original columns in the output of the liftover
let additional_columns = site_fields.join("\t");

// Return the formatted output string
return Some(format!(
Expand Down
10 changes: 5 additions & 5 deletions src/parse_gtf.rs
Original file line number Diff line number Diff line change
Expand Up @@ -368,22 +368,22 @@ pub fn read_annotation_file(file_path: &str, is_gtf: bool) -> Result<HashMap<Str
if feature.end < cds_start {
// 5'UTR for positive strand
transcript.utr5_len = Some(transcript.utr5_len.unwrap_or(0) + feature.length);
println!("5'UTR detected with length {}", feature.length);
// println!("5'UTR detected with length {}", feature.length);
} else if feature.start > cds_end {
// 3'UTR for positive strand
transcript.utr3_len = Some(transcript.utr3_len.unwrap_or(0) + feature.length);
println!("3'UTR detected with length {}", feature.length);
// println!("3'UTR detected with length {}", feature.length);
}
},
"-" => {
if feature.start > cds_end {
// 5'UTR for negative strand (logic is reversed)
transcript.utr5_len = Some(transcript.utr5_len.unwrap_or(0) + feature.length);
println!("5'UTR (neg strand) detected with length {}", feature.length);
// println!("5'UTR (neg strand) detected with length {}", feature.length);
} else if feature.end < cds_start {
// 3'UTR for negative strand (logic is reversed)
transcript.utr3_len = Some(transcript.utr3_len.unwrap_or(0) + feature.length);
println!("3'UTR (neg strand) detected with length {}", feature.length);
// println!("3'UTR (neg strand) detected with length {}", feature.length);
}
},
_ => {}
Expand Down Expand Up @@ -487,7 +487,7 @@ mod tests {
init(); // Initialize logging if necessary

// file path
let gtf_file_path = "/home/150/as7425/R2Dtool/test/gencode_v38.gtf";
let gtf_file_path = "./test/gencode_v38.gtf";

// Read the GTF file
let transcripts = read_annotation_file(gtf_file_path, true).expect("Failed to read GTF file");
Expand Down

0 comments on commit 1c4e079

Please sign in to comment.