Pipeline for ingesting high-coverage PB/ONT aligned BAMs and format them for re-processing #418
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,73 @@ | ||
| FROM python:3.9.16-slim-bullseye | ||
|
There was a problem hiding this comment. These dockerfiles are so long, and have so much repeated content. |
||
|
|
||
| COPY remove_duplicate_ont_aln.py /opt/ | ||
| COPY remove_duplicate_ont_namesorted_unaligned.py /opt/ | ||
|
|
||
| RUN pip install pysam==0.21.0 | ||
|
|
||
| # install gcloud and gsutil cli | ||
| ARG DEBIAN_FRONTEND=noninteractive | ||
| RUN apt-get -qqy update --fix-missing && \ | ||
| apt-get -qqy dist-upgrade && \ | ||
| apt-get -qqy install --no-install-recommends \ | ||
| apt-transport-https \ | ||
| ca-certificates \ | ||
| gnupg \ | ||
| zlib1g-dev \ | ||
| curl \ | ||
| wget \ | ||
| tree \ | ||
| tabix && \ | ||
| echo "deb [signed-by=/usr/share/keyrings/cloud.google.gpg] https://packages.cloud.google.com/apt cloud-sdk main" | tee -a /etc/apt/sources.list.d/google-cloud-sdk.list && \ | ||
| curl https://packages.cloud.google.com/apt/doc/apt-key.gpg | apt-key --keyring /usr/share/keyrings/cloud.google.gpg add - && \ | ||
| apt-get -qqy update && \ | ||
| apt-get -qqy install --no-install-recommends google-cloud-cli && \ | ||
| gcloud config set core/disable_usage_reporting true && \ | ||
| gcloud config set component_manager/disable_update_check true && \ | ||
| gcloud config set metrics/environment github_docker_image && \ | ||
| apt-get -qqy purge gnupg && \ | ||
| apt-get -qqy clean && \ | ||
| rm -rf /tmp/* \ | ||
| /var/tmp/* \ | ||
| /var/cache/apt/* \ | ||
| /var/lib/apt/lists/* \ | ||
| /usr/share/man/?? \ | ||
| /usr/share/man/??_* | ||
|
|
||
| # install latest samtools | ||
| ARG DEBIAN_FRONTEND=noninteractive | ||
| ARG SAMTOOLS_VERSION=1.18 | ||
| ARG BCFTOOLS_VERSION=1.18 | ||
| RUN apt-get -qqy update --fix-missing && \ | ||
| apt-get -qqy dist-upgrade && \ | ||
| apt-get -qqy install --no-install-recommends \ | ||
| ca-certificates \ | ||
| libbz2-dev \ | ||
| libcurl4-openssl-dev \ | ||
| liblzma-dev \ | ||
| libncurses5-dev \ | ||
| autoconf \ | ||
|
There was a problem hiding this comment. Note that all this stuff -- automake, gcc, wget, etc. -- is left in the final image. There was a problem hiding this comment. I now see that you remove most of it below. I still think it would be cleaner to remove it via staging. |
||
| automake \ | ||
| bzip2 \ | ||
| gcc \ | ||
| make \ | ||
| wget \ | ||
| zlib1g-dev && \ | ||
| wget https://github.com/samtools/samtools/releases/download/${SAMTOOLS_VERSION}/samtools-${SAMTOOLS_VERSION}.tar.bz2 && \ | ||
| tar xjf samtools-${SAMTOOLS_VERSION}.tar.bz2 && \ | ||
| cd samtools-${SAMTOOLS_VERSION} && ./configure --without-curses --enable-libcurl && make -s all all-htslib && make install install-htslib && cd - && \ | ||
| rm -rf samtools-${SAMTOOLS_VERSION}* && \ | ||
| wget https://github.com/samtools/bcftools/releases/download/${BCFTOOLS_VERSION}/bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ | ||
| tar xjf bcftools-${BCFTOOLS_VERSION}.tar.bz2 && \ | ||
| cd bcftools-${BCFTOOLS_VERSION} && ./configure --without-curses && make -s && make install && cd - && \ | ||
| rm -rf bcftools-${BCFTOOLS_VERSION}* && \ | ||
| apt-get -qqy purge autoconf automake bzip2 gcc make wget && \ | ||
| apt-get -qqy clean && \ | ||
| rm -rf /tmp/* \ | ||
| /var/tmp/* \ | ||
| /var/cache/apt/* \ | ||
| /var/lib/apt/lists/* \ | ||
| /usr/share/man/?? \ | ||
| /usr/share/man/??_* && \ | ||
| samtools --help && \ | ||
| bcftools --help | ||
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,12 @@ | ||
| VERSION = 0.1.2 | ||
| TAG1 = us.gcr.io/broad-dsp-lrma/lr-bam-dedup:$(VERSION) | ||
| TAG2 = us.gcr.io/broad-dsp-lrma/lr-bam-dedup:latest | ||
|
|
||
| all: build push | ||
|
|
||
| build: | ||
| docker build -t $(TAG1) -t $(TAG2) . | ||
|
|
||
| push: | ||
| docker push $(TAG1) | ||
| docker push $(TAG2) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,76 @@ | ||
| import argparse | ||
| import pysam | ||
|
|
||
|
|
||
| def main(): | ||
| parser = argparse.ArgumentParser(description='Remove redundant alignment records from ONT BAM file', | ||
| prog='remove_redundant_reads') | ||
| parser.add_argument('-p', '--prefix', type=str, default="shard", help="Output prefix") | ||
| parser.add_argument('-a', '--annotations', type=str, help="Annotations on (potential) duplicate reads") | ||
|
|
||
| parser.add_argument('bam', type=str, help="BAM") | ||
| args = parser.parse_args() | ||
|
|
||
| # create a dict of set's, a trick to avoid Hash collisions | ||
| guilty_dict_per_chr = dict() | ||
| with open(args.annotations) as f: | ||
| for line in f: | ||
| arr = line.strip().split('\t') | ||
| name = arr[0] | ||
| chrom = arr[2] | ||
| guilty_dict_per_chr.setdefault(chrom, set()) | ||
| guilty_dict_per_chr[chrom].add(name) | ||
|
|
||
| print("chromosomes on which there are duplicate records:") | ||
| print(f" {guilty_dict_per_chr}") | ||
|
|
||
| # Silence message about the .bai file not being found. | ||
| pysam.set_verbosity(0) | ||
|
|
||
| num_alignments, num_dropped_alignments = 0, 0 | ||
| duplicate_signatures = [] | ||
| bf = pysam.Samfile(args.bam, 'rb', check_sq=False) | ||
| with pysam.Samfile(f'{args.prefix}.bam', 'wb', header=bf.header) as out: | ||
| # we rely on the observation that for coordinate sorted BAM, | ||
| # duplicate records will appear in blocks, hence once we step off a position with duplicates, we start afresh | ||
| current_position = -1 | ||
| current_signatures = set() | ||
| for read in bf: | ||
| num_alignments += 1 | ||
|
|
||
| chrom = read.reference_name | ||
| n = read.query_name | ||
| if chrom in guilty_dict_per_chr and n in guilty_dict_per_chr[chrom]: | ||
|
|
||
| mq = read.mapping_quality | ||
| sam_flag = read.flag | ||
| pos = read.reference_start | ||
| cigar = read.cigarstring | ||
| signature = f"{n}-{chrom}-{pos}-{mq}-{sam_flag}-{cigar}" | ||
|
There was a problem hiding this comment. I note that you've added the (very lengthy!) cigar string to the signature. I believe that the signature is already unnecessarily lengthy. |
||
|
|
||
| if current_position != pos: # new position, let's write and reset | ||
| out.write(read) | ||
| current_position = pos | ||
| current_signatures = set() | ||
| current_signatures.add(signature) | ||
| elif signature in current_signatures: # You're a duplicate record, and not appearing for the 1st time! | ||
| num_dropped_alignments += 1 | ||
| duplicate_signatures.append(signature) # same signature may appear more than twice, hence list and append | ||
| pass | ||
| else: # you are in a new group of duplicates that map to this location | ||
| out.write(read) | ||
| current_signatures.add(signature) | ||
| else: | ||
| out.write(read) | ||
|
|
||
| print(f'num_alignments: {num_alignments}') | ||
| print(f'num_dropped_alignments: {num_dropped_alignments}') | ||
| print(f'num_kept_alignments: {num_alignments - num_dropped_alignments}') | ||
|
|
||
| with open(f'{args.prefix}.duplicate.signatures.txt', 'w') as out: | ||
| for sig in duplicate_signatures: | ||
| out.write(f"{sig}\n") | ||
|
|
||
|
|
||
| if __name__ == "__main__": | ||
| main() | ||
Oops, something went wrong.
Oops, something went wrong.
Add this suggestion to a batch that can be applied as a single commit.
This suggestion is invalid because no changes were made to the code.
Suggestions cannot be applied while the pull request is closed.
Suggestions cannot be applied while viewing a subset of changes.
Only one suggestion per line can be applied in a batch.
Add this suggestion to a batch that can be applied as a single commit.
Applying suggestions on deleted lines is not supported.
You must change the existing code in this line in order to create a valid suggestion.
Outdated suggestions cannot be applied.
This suggestion has been applied or marked resolved.
Suggestions cannot be applied from pending reviews.
Suggestions cannot be applied on multi-line comments.
Suggestions cannot be applied while the pull request is queued to merge.
Suggestion cannot be applied right now. Please check back later.
There was a problem hiding this comment.
Choose a reason for hiding this comment
The reason will be displayed to describe this comment to others. Learn more.
I haven't looked at this in detail. I'm going to assume it's OK.