Preprocessing some human samples.

1 Introduction

export HS_TYPE=gene
export HS_TAG=ID

1.1 Testing directory

I think for speed sake, I will run these through salmon.

I have been making some changes to my pipeline recently; so I will put 100k reads into a test directory.

cd preprocessing
mkdir test
cd test
less ../SRR/*R1* | head -n 400000 > r1.fastq
less ../SRR/*R2* | head -n 400000 > r2.fastq
gzip *.fastq
inputs=$(/bin/ls outputs/12fastp/*-fastp.fastq.xz | tr '\n' ':')
cyoa --method salmon --species hg38_111 --input $inputs --libtype CDS --jprefix 20

1.2 Trimming via fastp

I am going to use that tree to run everything first, notably I want my tools to collect more thorough statistics on runtime etc to improve my nascent heuristics to choose memory/time on the cluster.

start=$(pwd)
for i in $(/bin/ls -d SR*); do
    cd ${start}/${i}
    mkdir unprocessed
    mv *.fastq.gz unprocessed/
    cyoa --method fastp --input $(/bin/ls -d unprocessed/*.fastq.gz | tr '\n' ':')
done
cd $start

1.3 Salmon quantification

start=$(pwd)
for i in $(/bin/ls -d SR*); do
    cd ${start}/${i}
    inputs=$(/bin/ls outputs/12fastp/*-fastp.fastq.xz | tr '\n' ':')
    cyoa --method salmon --species hg38_111 --input $inputs --libtype CDS --jprefix 20
done
cd $start

2 Load annotations

hs_annot <- load_biomart_annotations(year = 2021, month = 02)

## The biomart annotations file already exists, loading from it.

annot <- hs_annot[["gene_annotations"]]

3 Collect preprocessing information

meta <- gather_preprocessing_metadata("sample_sheets/PRJNA675090.xlsx")

## Did not find the condition column in the sample sheet.

## Filling it in as undefined.

## Did not find the batch column in the sample sheet.

## Filling it in as undefined.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## Writing new metadata to: sample_sheets/PRJNA675090_modified.xlsx

## Deleting the file sample_sheets/PRJNA675090_modified.xlsx before writing the tables.

4 Create SE

tx_gene_map <- hs_annot[["gene_tx_map"]]
tx_gene_map[["transcript"]] <- gsub(x = tx_gene_map[["transcript"]], pattern = "\\.\\d+$", replacement = "")
rownames(tx_gene_map) <- make.names(gsub(x = rownames(tx_gene_map), pattern = "\\.\\d+$", replacement = ""), unique = TRUE)

## Error in `rownames<-`(`*tmp*`, value = character(0)): attempt to set 'rownames' on an object with no dimensions

hs_se_tx <- create_se(meta[["new_meta"]], gene_info = annot, file_column = "salmon_count_table")

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 41 rows(samples) and 22 columns(metadata fields).

## Warning in create_se(meta[["new_meta"]], gene_info = annot, file_column = "salmon_count_table"): Some
## samples were removed when cross referencing the samples against the count data.

## Warning in create_se(meta[["new_meta"]], gene_info = annot, file_column = "salmon_count_table"): Even
## after changing the rownames in gene info, they do not match the count table.

## Even after changing the rownames in gene info, they do not match the count table.

## Here are the first few rownames from the count tables:

## ENST00000000233.10, ENST00000000412.8, ENST00000000442.11, ENST00000001008.6, ENST00000001146.7, ENST00000002125.9

## Here are the first few rownames from the gene information table:

## ENSG00000000003, ENSG00000000005, ENSG00000000419, ENSG00000000457, ENSG00000000460, ENSG00000000938

## Bringing together the count matrix and gene information.

## Some annotations were lost in merging, setting them to 'undefined'.

## Saving the summarized experiment to 'se.rda'.

## The final summarized experiment has 97117 rows and 22 columns.

hs_se_gene <- create_se(meta[["new_meta"]], gene_info = annot, file_column = "salmon_count_table",
                        tx_gene_map = tx_gene_map)

## Reading the sample metadata.

## Checking the state of the condition column.

## Checking the state of the batch column.

## Checking the condition factor.

## The sample definitions comprises: 41 rows(samples) and 22 columns(metadata fields).

## In some cases, (notably salmon) the format of the IDs used by this can be tricky.
## It is likely to require the transcript ID followed by a '.' and the ensembl column:
## 'transcript_version', which is explicitly different than the gene version column.
## If this is not correctly performed, very few genes will be observed

## reading in files with read_tsv

## 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

## Error in `colnames<-`(`*tmp*`, value = c("tx", "gene")): attempt to set 'colnames' on an object with less than two dimensions

5 Set condition/batch

I have extracted two potentially interesting columns from the metadata, after that I will need to read more carefully in the paper to try to get a sense of what is what…

hs_se <- set_se_conditions(hs_se_gene, fact = "controlp") %>%
  set_se_batches(fact = "type")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colData': error in evaluating the argument 'x' in selecting a method for function 'colData': object 'hs_se_gene' not found

6 A couple plots

plot_legend(hs_se)

## Error: object 'hs_se' not found

plot_libsize(hs_se)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_libsize': object 'hs_se' not found

plot_nonzero(hs_se)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_nonzero': object 'hs_se' not found

Presumably we will need to separate the exome and rnaseq data; but for the moment I will leave them together to see what I can see.

hs_norm <- normalize(hs_se, transform = "log2", convert = "cpm", filter = TRUE, norm = "quant")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'input' in selecting a method for function 'normalize': object 'hs_se' not found

plot_pca(hs_norm)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_pca': object 'hs_norm' not found

hs_nb <- normalize(hs_se,  transform = "log2", convert = "cpm", filter = TRUE, batch = "combat")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'input' in selecting a method for function 'normalize': object 'hs_se' not found

plot_pca(hs_nb)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_pca': object 'hs_nb' not found

Ok, nevermind, the exome data is too different.

hs_rna <- subset_se(hs_se, subset = "batch=='rnaseq'")

## Error: object 'hs_se' not found

hs_rna_norm <- normalize(hs_rna, transform = "log2", convert = "cpm", filter = TRUE, norm = "quant")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'input' in selecting a method for function 'normalize': object 'hs_rna' not found

norm_pca <- plot_pca(hs_rna_norm)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_pca': object 'hs_rna_norm' not found

pp(file = "images/norm_pca.pdf")
norm_pca

## Error: object 'norm_pca' not found

dev.off()

## png 
##   2

##hs_rna_nb <- normalize(hs_rna, transform = "log2", convert = "cpm", filter = TRUE, batch = "svaseq")
##plot_pca(hs_rna_nb)

7 Exclude odd sample

In our initial PCA plot, SRR12999746 looks particularly strange.

hs_idx <- colnames(hs_rna) != "SRR12999746"

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colnames': object 'hs_rna' not found

kept_ids <- colnames(hs_rna)[hs_idx]

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'colnames': object 'hs_rna' not found

hs_excluded <- subset_se(hs_rna, ids = kept_ids)

## Error: object 'hs_rna' not found

excluded_norm <- normalize(hs_excluded, filter = TRUE, convert = "cpm", transform = "log2")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'input' in selecting a method for function 'normalize': object 'hs_excluded' not found

plot_pca(excluded_norm)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'plot_pca': object 'excluded_norm' not found

8 First try DE

hs_de <- all_pairwise(hs_excluded, filter = TRUE, model_svs = "svaseq",
                      model_fstring = "~ 0 + condition", force = TRUE)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'object' in selecting a method for function 'pData': object 'hs_excluded' not found

hs_table <- combine_de_tables(hs_de, excel = "excel/control_vs_hcm.xlsx")

## Deleting the file excel/control_vs_hcm.xlsx before writing the tables.

## Error in h(simpleError(msg, call)): error in evaluating the argument 'expt' in selecting a method for function 'colors': object 'hs_de' not found

hs_sig <- extract_significant_genes(hs_table, excel = "excel/control_vs_hcm-sig.xlsx")

## Deleting the file excel/control_vs_hcm-sig.xlsx before writing the tables.

## Error: object 'hs_table' not found

9 Ontology shenanigans

hs_gp <- all_cprofiler(hs_sig, hs_table)

## Error: object 'hs_sig' not found

A few specific mSigDB categories: C8, C7, C5, C3, C2

This is actually C2 because I need to download a new copy of msigdb in the new format.

9.1 New copy of msigdb

I downloaded a fresh copy of the human mSigDB because my copy was in the old format and that no longer parses well.

9.1.1 Load msigdb

up_cp_c8 <- simple_clusterprofiler(input_up, de_table = table,
                                   msig_db = "reference/msigdb_v2024.1.Hs.db",
                                   do_msigdb = TRUE, msigdb_category = "C8",
                                   do_mesh = TRUE, do_dose = TRUE, orgdb_from = "ENSEMBL")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'input_up' not found

9.2 C2

c2 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c2",
                          id_type = "entrez")
c2_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c2)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'sig_genes' in selecting a method for function 'simple_cp_enricher': object 'input_up' not found

written <- write_xlsx(data = as.data.frame(c2_enricher), excel = "excel/c2_result.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'c2_enricher' not found

9.3 C3

c3 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c3",
                          id_type = "entrez")
c3_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c3)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'sig_genes' in selecting a method for function 'simple_cp_enricher': object 'input_up' not found

written <- write_xlsx(data = as.data.frame(c3_enricher), excel = "excel/c3_result.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'c3_enricher' not found

9.4 C4

c4 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c4",
                          id_type = "entrez")
c4_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c4)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'sig_genes' in selecting a method for function 'simple_cp_enricher': object 'input_up' not found

written <- write_xlsx(data = as.data.frame(c4_enricher), excel = "excel/c4_result.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'c4_enricher' not found

9.5 C5

c5 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c5",
                          id_type = "entrez")
c5_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c5)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'sig_genes' in selecting a method for function 'simple_cp_enricher': object 'input_up' not found

written <- write_xlsx(data = as.data.frame(c5_enricher), excel = "excel/c5_result.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'c5_enricher' not found

9.6 C6

C6 returns no hits.

c6 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c6",
                          id_type = "entrez")
c6_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c6)
written <- write_xlsx(data = as.data.frame(c6_enricher), excel = "excel/c6_result.xlsx")

9.7 C7

c7 <- load_gmt_signatures(signatures = "reference/msigdb_v2024.1.Hs.db", signature_category = "c7",
                          id_type = "entrez")
c7_enricher <- simple_cp_enricher(sig_genes = input_up, de_table = table, db = c7)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'sig_genes' in selecting a method for function 'simple_cp_enricher': object 'input_up' not found

written <- write_xlsx(data = as.data.frame(c7_enricher), excel = "excel/c7_result.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'c7_enricher' not found

input_up <- as.data.frame(hs_sig[["deseq"]][["ups"]][[1]])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'hs_sig' not found

input_down <- as.data.frame(hs_sig[["deseq"]][["downs"]][[1]])

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'hs_sig' not found

table <- hs_table[["data"]][[1]]

## Error: object 'hs_table' not found

up_gp <- simple_gprofiler(input_up)

## Error in h(simpleError(msg, call)): error in evaluating the argument 'table' in selecting a method for function '%in%': object 'input_up' not found

up_cp_c8 <- simple_clusterprofiler(input_up, de_table = table,
                                   do_msigdb = TRUE, msigdb_category = "C8",
                                   do_mesh = TRUE, do_dose = TRUE, orgdb_from = "ENSEMBL")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'x' in selecting a method for function 'rownames': object 'input_up' not found

write_xlsx(data = as.data.frame(up_cp_c8[["msigdb_data"]]), excel = "excel/c2_increased.xlsx")

## Error in h(simpleError(msg, call)): error in evaluating the argument 'data' in selecting a method for function 'write_xlsx': error in evaluating the argument 'x' in selecting a method for function 'as.data.frame': object 'up_cp_c8' not found