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 21 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 21 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 21 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 = "donor") %>%
  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

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

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

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

plot_pca(hs_rna_nb)

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