Skip to content

Commit

Permalink
Quick script to pull down HLA-A from a test BAM file
Browse files Browse the repository at this point in the history
  • Loading branch information
kvg committed Jun 13, 2024
1 parent 0c59876 commit a772e3b
Show file tree
Hide file tree
Showing 4 changed files with 24 additions and 10 deletions.
4 changes: 2 additions & 2 deletions Cargo.lock

Some generated files are not rendered by default. Learn more about how customized files appear on GitHub.

13 changes: 13 additions & 0 deletions go.sh
Original file line number Diff line number Diff line change
@@ -0,0 +1,13 @@
#!/bin/bash

set -euxo pipefail

cargo build --release >/dev/null 2>&1

./target/debug/hidive fetch \
-r \
-o HG002.HLA-A.bam \
-l chr6:29,940,532-29,947,870 \
gs://fc-bb12940d-f7ba-4515-8a98-42de84f63c34/HPRC_grch38/PBCCSWholeGenome/HG002/alignments/HG002.bam

samtools index HG002.HLA-A.bam
2 changes: 1 addition & 1 deletion src/skydive/src/stage.rs
Original file line number Diff line number Diff line change
Expand Up @@ -225,7 +225,7 @@ pub fn stage_data(
reads
.iter()
.for_each(|read| {
if !seen_read_ids.contains(read.qname()) && read_spans_locus(read.reference_start(), read.reference_end(), loci) {
if !seen_read_ids.contains(read.qname()) && (!require_spanning_reads || read_spans_locus(read.reference_start(), read.reference_end(), loci)) {
let _ = bam_writer.write(read);
}

Expand Down
15 changes: 8 additions & 7 deletions src/skydive/src/utils.rs
Original file line number Diff line number Diff line change
Expand Up @@ -29,28 +29,29 @@ pub fn parse_loci(loci_list: &Vec<String>) -> HashSet<(String, u64, u64)> {

pub fn parse_locus(locus: String) -> Result<(String, u64, u64)> {
let l_fmt = locus.replace(",", "");
let parts: Vec<&str> = l_fmt.split(|c| (c == ':' || c == '-')).collect();
let parts1: Vec<&str> = l_fmt.split(|c| c == ':').collect();
let parts2: Vec<&str> = parts1[1].split(|c| c == '-').collect();

let chr = parts[0].to_string();
let chr = parts1[0].to_string();

if parts.len() == 2 {
let start = match parts[1].parse::<u64>() {
if parts2.len() == 1 {
let start = match parts2[0].parse::<u64>() {
Ok(val) => val,
Err(e) => {
return Err(anyhow::Error::new(e));
}
};

Ok((chr, start - 1000, start + 1000))
} else if parts.len() == 3 {
let start = match parts[1].parse::<u64>() {
} else if parts2.len() == 2 {
let start = match parts2[0].parse::<u64>() {
Ok(val) => val,
Err(e) => {
return Err(anyhow::Error::new(e));
}
};

let stop = match parts[2].parse::<u64>() {
let stop = match parts2[1].parse::<u64>() {
Ok(val) => val,
Err(e) => {
return Err(anyhow::Error::new(e));
Expand Down

0 comments on commit a772e3b

Please sign in to comment.