finished star and cellranger

nf-core · Nov 10, 2022 · 1415f2e · 1415f2e
1 parent c4673d9
commit 1415f2e
Show file tree

Hide file tree

Showing 7 changed files with 41 additions and 10 deletions.
diff --git a/bin/cellranger_mtx_to_h5ad.py b/bin/cellranger_mtx_to_h5ad.py
@@ -3,7 +3,7 @@
 import argparse
 import os
 import pandas as pd
-import scipy
+from scipy import io
 from anndata import AnnData
 
 def mtx_to_adata(mtx_h5: str, sample: str, verbose: bool = False):
@@ -20,12 +20,24 @@ def mtx_to_adata(mtx_h5: str, sample: str, verbose: bool = False):
 
 def write_counts(
     adata: AnnData,
+    txp2gene: str,
     out: str,
     verbose: bool = True,):
 
-    pd.DataFrame(adata.var.index).to_csv(os.path.join(out, "features.tsv"), sep="\t", index=False, header=None)
+    features = pd.DataFrame()
+    features["id"] = adata.var.index
+
+    # if txp2gene file is available enrich features file with gene names
+    if txp2gene:
+        t2g = pd.read_table(f"{txp2gene}/star/geneInfo.tab", header=None, skiprows=1)
+        print(t2g)
+        id2name = {e[0]: e[1] for _, e in t2g.iterrows()}
+        print(id2name)
+        features["name"] = adata.var.index.map(id2name)
+
+    features.to_csv(os.path.join(out, "features.tsv"), sep="\t", index=False, header=None)
     pd.DataFrame(adata.obs.index).to_csv(os.path.join(out, "barcodes.tsv"), sep="\t", index=False, header=None)
-    scipy.io.mmwrite(os.path.join(out, "matrix.mtx"), adata.X.T, field="integer")
+    io.mmwrite(os.path.join(out, "matrix.mtx"), adata.X.T, field="integer")
 
     if verbose:
         print("Wrote features.tsv, barcodes.tsv, and matrix.mtx files to {}".format(args["out"]))
@@ -40,6 +52,7 @@ def write_counts(
     parser.add_argument("-v", "--verbose", dest="verbose", help="Toggle verbose messages", default=False)
     parser.add_argument("-s", "--sample", dest="sample", help="Sample name")
     parser.add_argument("-o", "--out", dest="out", help="Output path.")
+    parser.add_argument("--txp2gene", dest="txp2gene", help="Transcript to gene (t2g) file.", nargs='?', const='')
 
     args = vars(parser.parse_args())
 
@@ -54,6 +67,7 @@ def write_counts(
 
     write_counts(
         adata,
+        args["txp2gene"],
         args["sample"],
         verbose=args["verbose"]
     )

diff --git a/bin/mtx_to_h5ad.py b/bin/mtx_to_h5ad.py
@@ -3,7 +3,7 @@
 import pandas as pd
 import argparse
 import os
-import scipy
+from scipy import io
 from anndata import AnnData
 
 def mtx_to_adata(
@@ -33,6 +33,7 @@ def mtx_to_adata(
 def write_counts(
     adata: AnnData,
     txp2gene: str,
+    star_index: str,
     out: str,
     verbose: bool = True,):
 
@@ -48,9 +49,15 @@ def write_counts(
         id2name = {e[1]: e[2] for _, e in t2g.iterrows()}
         features["name"] = adata.var.index.map(id2name)
 
+    # if star_index file is available enrich features file with gene names
+    if star_index:
+        t2g = pd.read_table(f"{star_index}/geneInfo.tab", header=None, skiprows=1)
+        id2name = {e[0]: e[1] for _, e in t2g.iterrows()}
+        features["name"] = adata.var.index.map(id2name)
+
     features.to_csv(os.path.join(out, "features.tsv"), sep="\t", index=False, header=None)
     pd.DataFrame(adata.obs.index).to_csv(os.path.join(out, "barcodes.tsv"), sep="\t", index=False, header=None)
-    scipy.io.mmwrite(os.path.join(out, "matrix.mtx"), adata.X.T, field="integer")
+    io.mmwrite(os.path.join(out, "matrix.mtx"), adata.X.T, field="integer")
 
     if verbose:
         print("Wrote features.tsv, barcodes.tsv, and matrix.mtx files to {}".format(args["out"]))
@@ -68,6 +75,7 @@ def write_counts(
     parser.add_argument("-o", "--out", dest="out", help="Output path.")
     parser.add_argument("-a", "--aligner", dest="aligner", help="Which aligner has been used?")
     parser.add_argument("--txp2gene", dest="txp2gene", help="Transcript to gene (t2g) file.", nargs='?', const='')
+    parser.add_argument("--star_index", dest="star_index", help="Star index folder containing geneInfo.tab.", nargs='?', const='')
 
     args = vars(parser.parse_args())
 
@@ -90,6 +98,7 @@ def write_counts(
     write_counts(
         adata,
         args["txp2gene"],
+        args["star_index"],
         args["sample"],
         verbose=args["verbose"]
     )

diff --git a/modules/local/mtx_to_h5ad.nf b/modules/local/mtx_to_h5ad.nf
@@ -12,6 +12,7 @@ process MTX_TO_H5AD {
     // for each sample, the sub-folders and files come directly in array.
     tuple val(meta), path(inputs)
     path txp2gene
+    path star_index
 
     output:
     path "${meta.id}/*h5ad", emit: h5ad
@@ -45,6 +46,7 @@ process MTX_TO_H5AD {
     cellranger_mtx_to_h5ad.py \\
         --mtx filtered_feature_bc_matrix.h5 \\
         --sample ${meta.id} \\
+        --txp2gene ${txp2gene} \\
         --out ${meta.id}/${meta.id}_matrix.h5ad
     """
 
@@ -59,6 +61,7 @@ process MTX_TO_H5AD {
             --barcode *count/counts_unfiltered/\${input_type}.barcodes.txt \\
             --feature *count/counts_unfiltered/\${input_type}.genes.txt \\
             --txp2gene ${txp2gene} \\
+            --star_index ${star_index} \\
             --out ${meta.id}/${meta.id}_\${input_type}_matrix.h5ad ;
     done
     """
@@ -73,6 +76,7 @@ process MTX_TO_H5AD {
         --barcode $barcodes_tsv \\
         --feature $features_tsv \\
         --txp2gene ${txp2gene} \\
+        --star_index ${star_index} \\
         --out ${meta.id}/${meta.id}_matrix.h5ad
     """
 

diff --git a/subworkflows/local/align_cellranger.nf b/subworkflows/local/align_cellranger.nf
@@ -5,7 +5,6 @@
 include {CELLRANGER_MKGTF} from "../../modules/nf-core/cellranger/mkgtf/main.nf"
 include {CELLRANGER_MKREF} from "../../modules/nf-core/cellranger/mkref/main.nf"
 include {CELLRANGER_COUNT} from "../../modules/nf-core/cellranger/count/main.nf"
-include {MTX_TO_H5AD     } from "../../modules/local/mtx_to_h5ad.nf"
 
 // Define workflow to subset and index a genome region fasta file
 workflow CELLRANGER_ALIGN {
@@ -43,4 +42,5 @@ workflow CELLRANGER_ALIGN {
     emit:
         ch_versions
         cellranger_out  = CELLRANGER_COUNT.out.outs
+        txp2gene = cellranger_index
 }
diff --git a/subworkflows/local/mtx_conversion.nf b/subworkflows/local/mtx_conversion.nf
@@ -9,6 +9,7 @@ workflow MTX_CONVERSION {
     mtx_matrices
     samplesheet
     txp2gene
+    star_index
 
     main:
         ch_versions = Channel.empty()
@@ -18,7 +19,8 @@ workflow MTX_CONVERSION {
         //
         MTX_TO_H5AD (
             mtx_matrices,
-            txp2gene
+            txp2gene,
+            star_index
         )
 
         //

diff --git a/subworkflows/local/starsolo.nf b/subworkflows/local/starsolo.nf
@@ -51,6 +51,7 @@ workflow STARSOLO {
 
     emit:
     ch_versions
+    star_index  = star_index
     star_result = STAR_ALIGN.out.tab
     star_counts = STAR_ALIGN.out.counts
     for_multiqc = STAR_ALIGN.out.log_final.collect{it[1]}.ifEmpty([])

diff --git a/workflows/scrnaseq.nf b/workflows/scrnaseq.nf
@@ -168,6 +168,7 @@ workflow SCRNASEQ {
         )
         ch_versions = ch_versions.mix(STARSOLO.out.ch_versions)
         ch_mtx_matrices = ch_mtx_matrices.mix(STARSOLO.out.star_counts)
+        ch_star_index = STARSOLO.out.star_index
         ch_multiqc_star = STARSOLO.out.for_multiqc
     }
 
@@ -181,15 +182,15 @@ workflow SCRNASEQ {
         )
         ch_versions = ch_versions.mix(CELLRANGER_ALIGN.out.ch_versions)
         ch_mtx_matrices = ch_mtx_matrices.mix(CELLRANGER_ALIGN.out.cellranger_out)
+        ch_txp2gene = CELLRANGER_ALIGN.out.txp2gene
     }
 
     // Run mtx to h5ad conversion subworkflow
     MTX_CONVERSION (
         ch_mtx_matrices,
         ch_input,
-        ch_txp2gene
-        // ch_genome_fasta,
-        // ch_filter_gtf
+        ch_txp2gene,
+        ch_star_index
     )
 
     //Add Versions from MTX Conversion workflow too