Metadata Sheet Discussion & Multilane merging

[apeltzer]

I'd like to have an optional way to input a custom design via a TSV formatted metadata sheet. The idea is, that this will be used in all cases. When you don't supply this TSV file in the future, it will be created based on your regular --reads input, of course, this means with missing values (thus skipping e.g. downstream popgen tools analysis).

Idea for metadata, basically based on what Sarek does:

ID  LabID   SeqType Organism    Treatment    R1  R2  Population (ISO3)   Group   
ABRAKADABRA JK1945  paired  UDG+ path/to/R1  path/to/R2  GER LBK
ABRAKADABRA JK1946  single  UDG- path/to/R1      GER LBK
HOKUSPOKUS  JK1948  paired  UDG0.5 path/to/R1  path/to/R2  TRY Anatolijan_Farmer

Comments ?

[maxulysse]

here's Sarek docs about our the TSV file input metadata: https://github.com/SciLifeLab/Sarek/blob/master/docs/INPUT.md

[apeltzer]

Okay, lets add BAM and BAI as well so we can start from that as well 👍

[jfy133]

From my experience, I think we need to have 'grouping' columns for:

1. Grouping per read pair (merging read pairs together)
2. Grouping per lane (concatenating lane-wise after adaptor removal, post-AR FASTQC running prior concatenation)
3. Grouping per sample (concatenating mapped BAMs of the same sample that has multiple libraries)

Does that make sense?

[jfy133]

An alternative would be to replace point 2 with sequencer - and we concatenate post-AR FASTQs only if the grouped sample is assigned to e.g. NextSeq (where you have a single injection onto the flowcell), and if assigned to HiSeq then merging of BAMs happen after mapping.

[apeltzer]

All three points are possible.

If you want to have 1., simply add a different ID for each line - then each line is treat indepdently :-)
If you want to have 2, simply add the same ID for each line and multiple internal samples. These will then be merged together after mapping individually to have all indiivdual lanes in one sample :-)
Is also possible via a join - Sarek does that as well :-)

[apeltzer]

ID  LabID   SeqType Organism    Treatment    R1  R2  Population (ISO3)   Group   BAM    BAI 
ABRAKADABRA JK1945  paired Human UDG+ path/to/R1  path/to/R2  GER LBK path/to/bam path/to/bai
ABRAKADABRA JK1946  single Human UDG- path/to/R1      GER LBK path/to/bam path/to/bai
HOKUSPOKUS  JK1948  paired Human UDG0.5 path/to/R1  path/to/R2  TRY Anatolijan_Farmer path/to/bam path/to/bai

ID = Identifier, basically the sample identifier to which multiple LabIDs might belong to. Samples (see JK1945, JK1946) are for example a paired end sequencing sample and a single end sequencing sample that have been sequenced from the same bone sample "ABRAKADABRA" and can be merged together. If users want to keep this separate, they can specify different IDs - if they use the same ID, the samples are merged together after initial mapping etc pp.

LabID = Internal sample id, e.g. sequencing library id. Multiple LabIDs can belong to the same "bone" or similar sample.

SeqType = Paired or Single end

Organism = might help in choosing a genome once people use iGenomes or similar resources at some point. Might not be required atm.

Treatment = UDG treatment Type

R1 = Forward read path
R2 = reverse Read path

Populatino ISO 3 = Population ISO3 Code, to choose for population genetics analysis (optional).
Group = Arbitrary free text grouping information, can be used to group a set of samples in pop gen (optional)
BAM = Path to BAM file, will ignore R1/R2 then and use this for popgen or downstream analysis from scratch (remap the file)
BAI = index required for doing the same

[jfy133]

What about this?

Sample_Name	Library_ID	Lane	SeqType	Sequencer	Organism	Strandedness	UDG_Treatment	R1	R2	BAM	BAI	Custom_Col_1	Custom_Col_2
Sample_1	Sample_1_Lib_1	1	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_1	2	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_1	3	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_1	4	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_2	1	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_2	2	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_2	3	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_1	Sample_1_Lib_2	4	pairedend	nextseq	NA	double	full	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_2	Sample_2_Lib_1	1	pairedend	hiseq	NA	double	half	/path/to/	/path/to/	/path/to/	/path/to/	Population_1	Group_1
Sample_2	Sample_2_Lib_2	2	singleend	hiseq	NA	double	half	/path/to/	NA	/path/to/	/path/to/	Population_1	Group_1

Where each row of the table is adapter removal'ed together (i.e. anything with the same Sample_ID, Library ID and Lane).

Then:

• If the sequencer is NextSeq (/MiSeq?), then all FASTQs after adapter removal are concatenated together.
• If the sequencer is HiSeq then all FASTQS are mapped independently, and mapped BAMs are merged (and unmapped bams also merged together independently).

The justification of the splitting by machine is that a given library that is split across multiple lanes comes from a single injection (i.e. single pipetting from the library aliquot) and represents a single library 'sampling', whereas for HiSeq, each 'lane' is from a different injection (i.e. different pipetting of the library aliquot) thus represents different 'samplings' from the library aliquot.

However all final (collapsed/concatenated/merged) files and MultiQC report info's are named by the 'Sample_Name'.

Custom columns are defined as you like (e.g. the population/ISO3 stuff), which EAGER2 pick ups if they are required for specific modules

[apeltzer]

I don't see the usefulness for sequencer but other than that I'm fine with the proposal 👍

Custom Fields are fine too - but I'll have to define some that can be used by users and not allow arbitrary ones ;-)

[drpatelh]

As you already know I am a big advocate of design files as input because they are able to capture all of the information required for the pipeline. 😉 However, I think we should also maybe have a general discussion as to how to define the fields in tsv files. This is something that could become quite disparate across pipelines, and it would be good to try and keep the convention for standard fields/headings the same across as many pipelines as possible.

Heres one I am planning for the chipseq pipeline:
https://github.com/drpatelh/chipseq/blob/master/docs/usage.md#--design

and the one I use for the atacseq pipeline:
https://github.com/nf-core/atacseq/blob/master/docs/usage.md#--design

This works if you sequence the same library more than once because you can merge based on the replicate column.

Also, I was wondering whether its absolutely necessary to have a .bai column because these normally reside with the .bam files anyway? Im assuming you are trying to account for this not being the case?

[apeltzer]

I agree that having this kind of streamlined is a good idea. The ID/group information is fine to be used, but the group here would be for adding some downstream metadata, e.g. which population group a sample belonged to. So we'd need that as an extra column and not like it's used in chipseq/atacseq. Still, the approach would be very similar and should work 👍

With respect to having lane IDs, I don't think that's necessary - if people want to keep these separately they can do that by having different IDs in the front.

[jfy133]

@apeltzer not necessarily. I realise I wasn't clear before - the lane library/lane/sequencer information is important for DeDup.

Lets say one is doing deep sequencing of a sample. You have multiple extracts and libraries from the same sample, one of these is sequenced on a NextSeq, another is sequenced on one lane of a HiSeq, and another on a different lane. As the DNA is from a single sample/individual you want the final output to be named with that (as indicated in ID in your original post).

However, to perform Dedup correctly (as intended i.e. as a PCR duplicate removal) - this should only be performed on a single library, not on different libraries/extracts of the same sample - because then this is not actually a PCR duplicate but a biological one.

Therefore, you need to indicate which libraries are OK to concatenate after adapterRemoval (with NextSeq), then you need to indicate which libraries (in this case lanes or Lab_ID) to perform DeDup on. but how this is performed depends on what type of 'lane' your flowcell has.

Obviously one could argue if it's still duplicated it's not 'useful' but it can be useful if the non-PCR related duplicate could be derived from two related species (e.g. in metagenomics), which you may want to utilise for abundance calculations etc.

[drpatelh]

Thats a valid point @jfy133 . However, if you want to make the most of this information you would need to populate the relevant fields in the Read Group for that sample. The simplest solution is to map each "run" of the sample individually adding the correct Read Group, merge libraries from the same sample, and then to re-perform duplicate marking and removal:
https://software.broadinstitute.org/gatk/documentation/article?id=6472

This will require quite a bit more to be implemented because the current set-up takes --reads as input which I assume means that multiple fastq files from the same sample are merged.

I have a similar implementation for the atacseq pipeline but I took a more crude approach where multiple runs of the same sample are given a unique read group ID, and subsequently merged.

[jfy133]

I see. I guess maybe Alex can take over discussion with the technical details here - if I understand it aDNA work does a slightly unusual series of steps compared to modern DNA because we often have such low coverage (but I'm an archaeologist who accidentally fell into bioinformatics, so what do I know ;) - thus I also may be getting confused in terminology at points).

I don't think the read group information is necessarily that important for the standard aDNA tools in EAGER (like in DeDup - correct me if I'm wrong @apeltzer). But if we are looking for a standardised 'design file' system then maybe that should be developed separately, then we 'hack' around it for aDNA contexts?

[apeltzer]

Hi both,

there are quite some different steps for aDNA, also with respect to duplicate removal which we typically don't do with MarkDuplicates but utilize DeDup (my own tool) that takes both start and end positions into account (read Peltzer et al 2016 for details on that). The problem is indeed, that DeDup cannot take read groups into account at the moment - so merging all based on read groups into a single BAM and then dedup based on read group information there isn't possible at the moment. The only prospect I see there is to dedup first and then merge the subsequent BAM files in all cases - which will be correct for both hiSeq and NextSeq output in my opinion.

I would like to be as similar as possible to other design files but if there are some specific parts missing for aDNA, we'd have to add these - chipseq/atacseq might have to at some point too 👍