The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.

Get started with biocohort

library(biocohort)
#> biocohort version 0.1.1
#> Cohort Objects for Subjects and Samples in Omics Studies
#> Documentation: https://www.samuelbharti.com/biocohort/

This article shows the path most studies take: a manifest file becomes a Cohort, the cohort answers questions about who has what sample, and it writes the sample sheet a pipeline expects.

From a manifest to a cohort

A manifest is one row per sample. Four columns carry the shape of the study: subject_id, assay, sample_id, and role. Every other column is metadata, either about the subject (species, genotype, sex, …) or about the sample (fastq_1, fastq_2, lane, …).

dir <- tempfile()
dir.create(dir)
writeLines(
  c(
    "subject_id,species,genotype,sex,assay,sample_id,role,fastq_1,fastq_2",
    "R1,rat,WT,F,wes,T1,tumor,t1_R1.fq.gz,t1_R2.fq.gz",
    "R1,rat,WT,F,wes,N1,normal,n1_R1.fq.gz,n1_R2.fq.gz",
    "R2,rat,KO,M,wes,T2,tumor,t2_R1.fq.gz,t2_R2.fq.gz",
    "R2,rat,KO,M,wes,N2,normal,n2_R1.fq.gz,n2_R2.fq.gz"
  ),
  file.path(dir, "manifest.csv")
)

parsed <- read_manifest(file.path(dir, "manifest.csv"))
parsed$subject_tbl
#> # A tibble: 2 × 4
#>   subject_id species genotype sex  
#>   <chr>      <chr>   <chr>    <chr>
#> 1 R1         rat     WT       F    
#> 2 R2         rat     KO       M
parsed$sample_map
#> # A tibble: 4 × 6
#>   subject_id assay sample_id role   fastq_1     fastq_2    
#>   <chr>      <chr> <chr>     <chr>  <chr>       <chr>      
#> 1 R1         wes   T1        tumor  t1_R1.fq.gz t1_R2.fq.gz
#> 2 R1         wes   N1        normal n1_R1.fq.gz n1_R2.fq.gz
#> 3 R2         wes   T2        tumor  t2_R1.fq.gz t2_R2.fq.gz
#> 4 R2         wes   N2        normal n2_R1.fq.gz n2_R2.fq.gz

read_manifest() reads every column as text, so an id like 007 keeps its leading zero. Build the cohort from the parsed tables:

cohort <- cohort_new(parsed$subject_tbl, parsed$sample_map)
cohort
#> 
#> ── Cohort
#> • 2 subjects (2 rat)
#> • 4 samples (4 wes)
#> ℹ Extra sample columns: fastq_1, fastq_2

A wide table instead

Some studies keep one row per subject, with one id column per assay. Turn that into the long form with manifest_from_wide() before building a cohort:

wide <- data.frame(
  subject_id = c("R1", "R2"),
  species = "rat",
  wes_tumor_id = c("T1", "T2"),
  wes_normal_id = c("N1", "N2"),
  stringsAsFactors = FALSE
)

id_cols <- data.frame(
  column = c("wes_tumor_id", "wes_normal_id"),
  assay = c("wes", "wes"),
  role = c("tumor", "normal"),
  stringsAsFactors = FALSE
)

manifest_from_wide(wide, id_cols)
#> # A tibble: 4 × 5
#>   subject_id species assay sample_id role  
#>   <chr>      <chr>   <chr> <chr>     <chr> 
#> 1 R1         rat     wes   T1        tumor 
#> 2 R2         rat     wes   T2        tumor 
#> 3 R1         rat     wes   N1        normal
#> 4 R2         rat     wes   N2        normal

Reading the cohort back

The accessors return plain tibbles, so the rest of a script can use ordinary dplyr code.

subjects(cohort)
#> # A tibble: 2 × 4
#>   subject_id species genotype sex  
#>   <chr>      <chr>   <chr>    <chr>
#> 1 R1         rat     WT       F    
#> 2 R2         rat     KO       M
samples(cohort, assay = "wes")
#> # A tibble: 4 × 6
#>   subject_id assay sample_id role   fastq_1     fastq_2    
#>   <chr>      <chr> <chr>     <chr>  <chr>       <chr>      
#> 1 R1         wes   T1        tumor  t1_R1.fq.gz t1_R2.fq.gz
#> 2 R1         wes   N1        normal n1_R1.fq.gz n1_R2.fq.gz
#> 3 R2         wes   T2        tumor  t2_R1.fq.gz t2_R2.fq.gz
#> 4 R2         wes   N2        normal n2_R1.fq.gz n2_R2.fq.gz
completeness(cohort, wide = TRUE)
#> # A tibble: 2 × 2
#>   subject_id   wes
#>   <chr>      <int>
#> 1 R1             2
#> 2 R2             2

cohort_filter() keeps a subset of subjects and returns a cohort that is still valid:

cohort_filter(cohort, genotype == "KO")
#> 
#> ── Cohort
#> • 1 subject (1 rat)
#> • 2 samples (2 wes)
#> ℹ Extra sample columns: fastq_1, fastq_2

Writing a pipeline sample sheet

sample_sheet() writes the sample list in the shape a pipeline expects. Built-in templates cover a few common nf-core pipelines:

sample_sheet_templates()
#> [1] "nf-core/rnaseq"  "nf-core/rnavar"  "nf-core/atacseq" "nf-core/sarek"
sample_sheet(cohort, template = "nf-core/sarek", assay = "wes")
#> # A tibble: 4 × 7
#>   patient sex   status sample  lane fastq_1     fastq_2    
#>   <chr>   <chr>  <int> <chr>  <int> <chr>       <chr>      
#> 1 R1      XX         1 T1         1 t1_R1.fq.gz t1_R2.fq.gz
#> 2 R1      XX         0 N1         1 n1_R1.fq.gz n1_R2.fq.gz
#> 3 R2      XY         1 T2         1 t2_R1.fq.gz t2_R2.fq.gz
#> 4 R2      XY         0 N2         1 n2_R1.fq.gz n2_R2.fq.gz

Registering an analysis

An AnalysisSpec records where an analysis writes its output and how to read it back. load_analysis() then resolves the path for every subject or pair and reads what it finds.

spec <- analysis_spec_new(
  name = "somatic_vars",
  assay = "wes",
  level = "pair",
  path_template = file.path(dir, "{pair_id}.tsv")
)
cohort <- analysis_register(cohort, spec)
analysis_list(cohort)
#> # A tibble: 1 × 6
#>   name         assay level format reader          root_key
#>   <chr>        <chr> <chr> <chr>  <chr>           <chr>   
#> 1 somatic_vars wes   pair  tsv    readr::read_tsv <NA>

Saving and loading a cohort

cohort_save(cohort, file.path(dir, "cohort.rds"))
reread <- cohort_read(file.path(dir, "cohort.rds"))
identical(subjects(cohort), subjects(reread))
#> [1] TRUE

A study with more than a manifest, a few paths, and a couple of analyses is easier to keep in one YAML file. See ?read_study_yaml for the file format.

Translating features across species

translate() moves a feature table from one species or genome build to another. Coordinate features go through a liftover backend. Gene features go through an ortholog backend. See ?translate and ?liftover_intervals for the full set of options, including how to register a custom backend.

Where to go next

These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.