Fix genome index generation; Reimplementation of count_matrix buildin…

…g; Fix column order matching

workflow/rules/coverage_tRNA.smk

@@ -18,7 +18,7 @@ rule mk_genome_tsv:
         mem_mb=2048,
     shell:
         """
-        samtools view -H {input} | grep SQ | cut -f 2,3 | sed -r 's/(SN|LN)://g' | tr " " "\t" > {output}
+        samtools view -H {input} | grep '@SQ' | cut -f 2,3 | sed -r 's/(SN|LN)://g' | tr " " "\t" > {output}
         """
 
 

workflow/scripts/build_tRNA_coverage_matrix_v1.R

@@ -16,13 +16,21 @@ print(coverage_files)
 coverage_files <- coverage_files[match(sample_sheet$name, names(coverage_files))]
 print(coverage_files)
 
-
 out <- coverage_files %>%
   map(read_tsv, col_names = FALSE) %>%
-  map(select, c(X4, X14)) %>%
-  purrr::reduce(inner_join, by = "X4") %>%
-  rename_with(function(x) c("Name", basename(as.character(coverage_files)))) 
-
-print(head(out), file = stderr())
-
-write_tsv(x = out, file = as.character(snakemake@output))
+  map(dplyr::select, c(X4, X14))
+
+ncol <- length(out)
+nrow <- nrow(out[[1]]) # assume that all elements in out have the same number of rows
+
+m <- matrix(0, nrow, ncol, dimnames = list(out[[1]]$X4, basename(as.character(coverage_files))))
+for (i in seq(length(out))) {
+    df <- out[[i]]
+    if (nrow(df) != nrow(m)) stop(sprintf("Incorrect dimensions! Element %d has %d rows instead of %d.",
+        i, nrow(df), nrow(m)))
+    #df <- df[match(rownames(m), df$X4),]
+    m[,i] <- df$X14
+}
+m <- as_tibble(m, rownames="Name")
+
+write_tsv(x = m, file = as.character(snakemake@output))

workflow/scripts/coverage_trna.sh

@@ -12,13 +12,13 @@ for C in $CHROMOSOMES; do
     awk -v j="$C" '$1~j' "${snakemake_input[annotation]}" >> $T
 done 
 
-bedtools sort -faidx "${snakemake_input[genome]}" -i $T > $SORTED 2> ${snakemake_log}
+bedtools sort -faidx "${snakemake_input[genome]}" -i $T > $SORTED
 
 bedtools coverage \
     -g "${snakemake_input[genome]}" \
     -sorted \
     -a $SORTED \
     -b "${snakemake_input[bam]}" | \
-    sort -k1,1 -k2,2n > "${snakemake_output}" 2>> ${snakemake_log}
+    sort -k1,1 -k2,2n > "${snakemake_output}"
 
 rm $SORTED

workflow/scripts/deseq2_trna_v1.R

@@ -27,8 +27,7 @@ sample_sheet <- read.csv(
     header=TRUE
 )
 
-
-colnames_order <- sapply(colnames(count_matrix), grep, x = sample_sheet$name )
+colnames_order <- match(sample_sheet$name,colnames(count_matrix))
 colnames(count_matrix)[colnames_order] <- rownames(sample_sheet)
 
 sample_sheet <- sample_sheet[match(colnames(count_matrix), rownames(sample_sheet)),]