diff --git a/src/bd_rhapsody/bd_rhapsody_sequence_analysis/test.py b/src/bd_rhapsody/bd_rhapsody_sequence_analysis/test.py index 8d7e68ec..57d699f0 100644 --- a/src/bd_rhapsody/bd_rhapsody_sequence_analysis/test.py +++ b/src/bd_rhapsody/bd_rhapsody_sequence_analysis/test.py @@ -32,6 +32,7 @@ reference_small_gtf = meta["resources_dir"] / "test_data" / "reference_small.gtf" reference_small_fa = meta["resources_dir"] / "test_data" / "reference_small.fa" bdabseq_panel_fa = meta["resources_dir"] / "test_data" / "BDAbSeq_ImmuneDiscoveryPanel.fasta" +sampletagsequences_fa = meta["resources_dir"] / "test_data" / "SampleTagSequences_HomoSapiens_ver1.fasta" config_file = Path("reference_config.yml") reference_file = Path("Rhap_reference.tar.gz") @@ -62,6 +63,8 @@ reference_fasta_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta")) with open(str(bdabseq_panel_fa), "r") as handle: bdabseq_panel_fasta_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta")) +with open(str(sampletagsequences_fa), "r") as handle: + sampletagsequences_fasta_dict = SeqIO.to_dict(SeqIO.parse(handle, "fasta")) # create in memory db reference_gtf_db = gffutils.create_db( @@ -306,6 +309,85 @@ def generate_bd_abc_fastq_files( return r1_reads, r2_reads +def generate_bd_smk_read( + cell_index: int = 0, + bead_version: str = "EnhV2", + umi_length: int = 14, + transcript_length: int = 72, + num_sample_tags: int = 3, +): + """ + Generate a BD Rhapsody SMK read pair for a given cell index. + + Args: + cell_index: The cell index to generate reads for. + bead_version: The bead version to use for generating the cell label. + umi_length: The length of the UMI to generate. + transcript_length: The length of the transcript to generate + num_sample_tags: The number of sample tags to use + + Returns: + A tuple of two strings, the first string being the R1 read and the second string being the R2 read. + """ + # generate metadata + per_row = np.floor((32967 - 1000) / 9) + per_col = np.floor((37059 - 1000) / 9) + + assert cell_index >= 0 and cell_index < per_row * per_col, f"cell_index must be between 0 and {per_row} * {per_col}" + x = 1000 + (cell_index % per_row) * 9 + y = 1000 + (cell_index // per_row) * 9 + instrument_id = "A00226" + run_id = "970" + flowcell_id = "H5FGVDMXY" + + meta_r1 = generate_bd_read_metadata(instrument_id=instrument_id, run_id=run_id, flowcell_id=flowcell_id, x=x, y=y, illumina_flag="1:N:0") + meta_r2 = generate_bd_read_metadata(instrument_id=instrument_id, run_id=run_id, flowcell_id=flowcell_id, x=x, y=y, illumina_flag="2:N:0") + + # generate r1 (cls1 + link + cls2 + link + cls3 + umi) + assert cell_index >= 0 and cell_index < 384 * 384 * 384 + cell_label = index_to_sequence(cell_index + 1, bead_version=bead_version) + # sample random umi + umi = "".join(random.choices("ACGT", k=umi_length)) + quality_r1 = "I" * (len(cell_label) + len(umi)) + r1 = f"{meta_r1}\n{cell_label}{umi}\n+\n{quality_r1}\n" + + # generate r2 by selecting the cell_index %% num_sample_tags sample tags + sampletag_index = cell_index % num_sample_tags + sampletag_seq = str(list(sampletagsequences_fasta_dict.values())[sampletag_index].seq) + smk_data = sampletag_seq[:transcript_length] + smk_suffix = "A" * (transcript_length - len(smk_data)) + quality_r2 = "I" * len(smk_data) + "#" * len(smk_suffix) + r2 = f"{meta_r2}\n{smk_data}{smk_suffix}\n+\n{quality_r2}\n" + + return r1, r2 + +def generate_bd_smk_fastq_files( + num_cells: int = 100, + num_reads_per_cell: int = 1000, + num_sample_tags: int = 3, +) -> Tuple[str, str]: + """ + Generate BD Rhapsody SMK FASTQ files for a given number of cells and transcripts per cell. + + Args: + num_cells: The number of cells to generate + num_reads_per_cell: The number of reads to generate per cell + num_sample_tags: The number of sample tags to use + + Returns: + A tuple of two strings, the first string being the R1 reads and the second string being the R2 reads. + """ + r1_reads = "" + r2_reads = "" + for cell_index in range(num_cells): + for _ in range(num_reads_per_cell): + r1, r2 = generate_bd_smk_read(cell_index, num_sample_tags=num_sample_tags) + r1_reads += r1 + r2_reads += r2 + + return r1_reads, r2_reads + +######################################################################################### # Prepare WTA, ABC, and SMK test data print("> Prepare WTA test data", flush=True) @@ -322,6 +404,13 @@ def generate_bd_abc_fastq_files( with gzip.open("ABCreads_R2.fq.gz", "wt") as f: f.write(abc_reads_r2_str) +print("> Prepare SMK test data", flush=True) +smk_reads_r1_str, smk_reads_r2_str = generate_bd_smk_fastq_files(num_cells=100, num_reads_per_cell=1000, num_sample_tags=3) +with gzip.open("SMKreads_R1.fq.gz", "wt") as f: + f.write(smk_reads_r1_str) +with gzip.open("SMKreads_R2.fq.gz", "wt") as f: + f.write(smk_reads_r2_str) + ######################################################################################### # Run executable @@ -330,14 +419,17 @@ def generate_bd_abc_fastq_files( subprocess.run([ meta['executable'], "--reads=WTAreads_R1.fq.gz;WTAreads_R2.fq.gz", - "--reads=ABCreads_R1.fq.gz;ABCreads_R2.fq.gz", f"--reference_archive={reference_file}", + "--reads=ABCreads_R1.fq.gz;ABCreads_R2.fq.gz", f"--abseq_reference={bdabseq_panel_fa}", + "--reads=SMKreads_R1.fq.gz;SMKreads_R2.fq.gz", + "--tag_names=1-Sample1;2-Sample2;3-Sample3", + "--sample_tags_version=human", "--output_dir=output", "--exact_cell_count=100", f"---cpus={meta['cpus'] or 1}", f"---memory={meta['memory_mb'] or 2048}mb", - "--output_seurat=seurat.rds", + # "--output_seurat=seurat.rds", "--output_mudata=mudata.h5mu", "--metrics_summary=metrics_summary.csv", "--pipeline_report=pipeline_report.html", @@ -351,11 +443,11 @@ def generate_bd_abc_fastq_files( assert (output_dir / "sample_Pipeline_Report.html").exists() assert (output_dir / "sample_RSEC_MolsPerCell_MEX.zip").exists() assert (output_dir / "sample_RSEC_MolsPerCell_Unfiltered_MEX.zip").exists() -assert (output_dir / "sample_Seurat.rds").exists() +# assert (output_dir / "sample_Seurat.rds").exists() assert (output_dir / "sample.h5mu").exists() # check individual outputs -assert Path("seurat.rds").exists() +# assert Path("seurat.rds").exists() assert Path("mudata.h5mu").exists() assert Path("metrics_summary.csv").exists() assert Path("pipeline_report.html").exists() @@ -365,17 +457,31 @@ def generate_bd_abc_fastq_files( assert data.n_obs == 100, "Number of cells is incorrect" assert "rna" in data.mod, "RNA data is missing" +assert "prot" in data.mod, "Protein data is missing" - +# check rna data data_rna = data.mod["rna"] assert data_rna.n_vars == 1, "Number of genes is incorrect" - -# row sum should be greater than 950 assert data_rna.X.sum(axis=1).min() > 950, "Number of reads per cell is incorrect" assert data_rna.var.Raw_Reads.sum() == 100000, "Number of reads is incorrect" -# TODO: add VDJ, SMK, ATAC, and targeted RNA to test - +# check prot data +data_prot = data.mod["prot"] +assert data_prot.n_vars == bdabseq_panel_fa.n_seqs, "Number of proteins is incorrect" +assert data_prot.X.sum(axis=1).min() > 950, "Number of reads per cell is incorrect" +assert data_prot.var.Raw_Reads.sum() == 100000, "Number of reads is incorrect" + +# check smk data +expected_sample_tags = (["SampleTag01_hs", "SampleTag02_hs", "SampleTag03_hs"] * 34)[:100] +expected_sample_names = (["Sample1", "Sample2", "Sample3"] * 34)[:100] +sample_tags = data_rna.obs["Sample_Tag"] +assert sample_tags.nunique() == 3, "Number of sample tags is incorrect" +assert sample_tags.tolist() == expected_sample_tags, "Sample tags are incorrect" +sample_names = data_rna.obs["Sample_Name"] +assert sample_names.nunique() == 3, "Number of sample names is incorrect" +assert sample_names.tolist() == expected_sample_names, "Sample names are incorrect" + +# TODO: add VDJ, ATAC, and targeted RNA to test ######################################################################################### diff --git a/src/bd_rhapsody/test_data/SampleTagSequences_HomoSapiens_ver1.fasta b/src/bd_rhapsody/test_data/SampleTagSequences_HomoSapiens_ver1.fasta new file mode 100644 index 00000000..3d5a42fa --- /dev/null +++ b/src/bd_rhapsody/test_data/SampleTagSequences_HomoSapiens_ver1.fasta @@ -0,0 +1,24 @@ +>SampleTag01_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGATTCAAGGGCAGCCGCGTCACGATTGGATACGACTGTTGGACCGG +>SampleTag02_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTGGATGGGATAAGTGCGTGATGGACCGAAGGGACCTCGTGGCCGG +>SampleTag03_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCGGCTCGTGCTGCGTCGTCTCAAGTCCAGAAACTCCGTGTATCCT +>SampleTag04_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGATTGGGAGGCTTTCGTACCGCTGCCGCCACCAGGTGATACCCGCT +>SampleTag05_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCTCCCTGGTGTTCAATACCCGATGTGGTGGGCAGAATGTGGCTGG +>SampleTag06_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTTACCCGCAGGAAGACGTATACCCCTCGTGCCAGGCGACCAATGC +>SampleTag07_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTGTCTACGTCGGACCGCAAGAAGTGAGTCAGAGGCTGCACGCTGT +>SampleTag08_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCCCCACCAGGTTGCTTTGTCGGACGAGCCCGCACAGCGCTAGGAT +>SampleTag09_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGTGATCCGCGCAGGCACACATACCGACTCAGATGGGTTGTCCAGG +>SampleTag10_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGCAGCCGGCGTCGTACGAGGCACAGCGGAGACTAGATGAGGCCCC +>SampleTag11_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCGCGTCCAATTTCCGAAGCCCCGCCCTAGGAGTTCCCCTGCGTGC +>SampleTag12_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGCCCATTCATTGCACCCGCCAGTGATCGACCCTAGTGGAGCTAAG diff --git a/src/bd_rhapsody/test_data/script.sh b/src/bd_rhapsody/test_data/script.sh index 19d97886..f8db0313 100644 --- a/src/bd_rhapsody/test_data/script.sh +++ b/src/bd_rhapsody/test_data/script.sh @@ -110,3 +110,32 @@ TTTGGAGGGTAGCTAGTTGCAGTTCGTGGTCGTTTC >IgD|IGHD|AHS0058|pAbO Catalog_940026 TGAGGGATGTATAGCGAGAATTGCGACCGTAGACTT EOF + +# this was obtained by running the command: +# docker run bdgenomics/rhapsody:2.2.1 cat /rhapsody/control_files/SampleTagSequences_HomoSapiens_ver1.fasta +cat > $OUT_DIR/SampleTagSequences_HomoSapiens_ver1.fasta <SampleTag01_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGATTCAAGGGCAGCCGCGTCACGATTGGATACGACTGTTGGACCGG +>SampleTag02_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTGGATGGGATAAGTGCGTGATGGACCGAAGGGACCTCGTGGCCGG +>SampleTag03_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCGGCTCGTGCTGCGTCGTCTCAAGTCCAGAAACTCCGTGTATCCT +>SampleTag04_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGATTGGGAGGCTTTCGTACCGCTGCCGCCACCAGGTGATACCCGCT +>SampleTag05_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCTCCCTGGTGTTCAATACCCGATGTGGTGGGCAGAATGTGGCTGG +>SampleTag06_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTTACCCGCAGGAAGACGTATACCCCTCGTGCCAGGCGACCAATGC +>SampleTag07_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGTGTCTACGTCGGACCGCAAGAAGTGAGTCAGAGGCTGCACGCTGT +>SampleTag08_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCCCCACCAGGTTGCTTTGTCGGACGAGCCCGCACAGCGCTAGGAT +>SampleTag09_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGTGATCCGCGCAGGCACACATACCGACTCAGATGGGTTGTCCAGG +>SampleTag10_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGCAGCCGGCGTCGTACGAGGCACAGCGGAGACTAGATGAGGCCCC +>SampleTag11_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGCGCGTCCAATTTCCGAAGCCCCGCCCTAGGAGTTCCCCTGCGTGC +>SampleTag12_hs|stAbO +GTTGTCAAGATGCTACCGTTCAGAGGCCCATTCATTGCACCCGCCAGTGATCGACCCTAGTGGAGCTAAG +EOF