Hey all. I am working a project to determine whether fibroblast cells still contain variants using whole genome sequencing (WGS) data from paired-end FASTQ files. my samplesheet_wgs.csv is: patient,sample,lane,fastq_1,fastq_2 p1,NORMAL,L001,p1_S_R1_001.fastq.gz,p1_S_R2_001.fastq.gz I am running nf-core/sarek v3.5.1 $ nextflow run nf-core/sarek -r 3.5.1 \ -profile singularity \ --input samplesheet_wgs.csv \ --outdir results_wgs \ --step mapping \ --genome GATK.GRCh38 \ --tools haplotypecaller,strelka,cnvkit,vep \ --vep_loftee \ -with-dag flowchart_wes.png I have two questions regarding how to complete the task: 1. Since the goal is to determine whether mutations are still present in the fibroblast sample, should I define the sample as "NORMAL" or "TUMOR" in the samplesheet_wgs.csv? 2. Because this is a WGS dataset, do I still need to include the following parameters, or are they optional? --vep_spliceai --spliceai_snv vep/spliceai_scores.raw.snv.hg38.vcf.gz \ --spliceai_snv_tbi vep/spliceai_scores.raw.snv.hg38.vcf.gz.tbi \ --spliceai_indel vep/spliceai_scores.raw.indel.hg38.vcf.gz \ --spliceai_indel vep/spliceai_scores.raw.indel.hg38.vcf.gz.tbi \ --vep_loftee \ Thank you very much for your guidance and help!

• Are you expecting somatic or germline mutations? If there is a chance of somatic variants, set the

status

column in the samplesheet to

1

(in that case also use

mutect2

instead of

haplotypecaller

). The default (i.e. no

status

column) is

0

, which uses germline methods. For Strelka to run without tumor-normal pairs, you have to use

0

, since it does not have a tumor-only mode (but you can use the germline results like tumor results, I think). •

--vep_spliceai

(and associated parameters) is definitely optional (not sure if there is any required parameter besides

--input

). Use it if you want to specifically flag variants that affect splicing. • You are running with just

-profile singularity

. Unless you have a single server (with lots of RAM), this will be highly inefficient. You need some configuration to tell Nextflow how to run the pipeline on your system (whatever that is -- cloud, HPC, etc.). See here for some tips.

Thank you very much for your guidance and support on how to run nf-core/sarek. My objective is to determine whether the fibroblast sample contains variants, regardless of whether they are of germline or somatic origin. To achieve this, I am considering running two separate analyses: 1. A germline variant calling workflow by setting the sample as “NORMAL” and using

haplotypecaller

. 2. A somatic variant calling workflow by setting the sample as “TUMOR” and using

mutect2

. Afterward, I would compare or combine the results from both analyses to determine whether variants are present in the fibroblast sample. Would this be an appropriate strategy for achieving my goal, or is there a more suitable approach? Thank you again for your support.

I have never done what you propose, but I you can definitely try that. Any sample will contain variants, it would be extremely unlikely that it conforms perfectly to the "reference genome". With no "normal control", you won't even know for sure a lot of times if a variant is germline or somatic. Because of that I would not be surprised if some of the variants show up with both callers. There is also Strelka and FreeBayes. FreeBayes has no built-in filter, so you have to do more work to remove artifacts (the author suggests using some "common-sense" filtering criteria). Are you looking for variants doing anything specific? If you find interesting candidates, definitely look at them in IGV (or something similar) to make sure they are not some weird mapping artifacts. Your WGS won't give you verry deep coverage, so it will be hard to find (somatic) variants with low allele frequency.

Thanks, Danhorn. My goal is to determine whether the fibroblast cells still contain any variants. Ideally, no variants should be detected, as these samples will be used as normal controls in subsequent analyses.

As I said, they will most likely contain germline variants (with respect to the reference genome). And sometimes it can be hard to distinguish somatic and germline variants, especially if the sequencing depth is not very high (because the sampling error would have a larger effect). 🤞

Thanks! I will plan to check both somatic and germline variants.

You are already doing this correctly with

--step variant_calling

. You just have to supply a samplesheet with CRAM (or BAM) files, not FASTQs. The ones to use are the "recalibrated" ones from the first run, and this samplesheet is automatically generated by the first run in

results/csv/recalibrated.csv

.

Thank you. First, I executed the following command to analyze the fibroblast lines WGS data sets: $ nextflow run nf-core/sarek -r 3.5.1 \ -profile singularity \ --input samplesheet_somatic.csv \ --outdir somatic \ --step mapping \ --genome GATK.GRCh38 \ --save_mapped \ --tools mutect2,cnvkit,vep \ --vep_loftee After that, I edited the file results/csv/recalibrated.csv, and modified line2 as follows: patient,sex,status,sample,cram,crai 1,N,0,tumor,tumor.recal.cram,tumor.recal.cram.crai I changed "1,N,1" to "1,N,0" because haplotypecaller can only run on sample labeled as normal. The I ran the following command:

nextflow run nf-core/sarek -r 3.5.1 \

-profile singularity \ --input recalibrated.csv \ --outdir germline \ --step variant_calling \ --genome GATK.GRCh38 \ --tools haplotypecaller,vep \ --vep_loftee The running time for "--step variant_calling" was approximately one-third of the time required for "--step mapping". For each sample, around 4M variants were generated from mutect2 and haplotypecaller. My questions is What methods or parameters in the Sarek pipeline can be used to filter and retain only high-confidence(true) variants, especially real somatic mutations?

Sarek itself does not do any filtering. How to best do that really depends on your exact question and the relative "cost" treating a variant as real when it is really an artifact, compared to discarding a variant that looks like an artifact but is real. Only you can answer these questions for your particular experiment. Some very basic things to do is to only use variants with a

FILTER

of

PASS

and to make sure you have enough read depth to make a reasonably reliable call. I don't have any real experience with germline variants, but for somatic variants you have to look at how well supported a variant is -- if there is a single mutated read, how do you know that's not a sequencing error? Can you see the variant on both strands? (Sometime there are too few reads, but if there are several on one strand and none on the other, I would suspect an artifact.) There are other criteria, and a lot of these are already built into the

FILTER

results. For somatic variants you can also filter by variant allele frequency (VAF) -- you will have to decide if what level you consider "actionable", i.e. if you are you interested in rare variants or not, or want to exclude very frequent ones, etc. I do most of these type of thing with

bcftools

. With VEP, you can use the

+split-vep

plugin to extract all of the annotation data and filter on that as well (i.e. you could remove silent mutations, or rank be severity, etc.). There is a new pipeline in development, #C05URH17X5Z, which you could try. I have not used it (and since it's in development, it might not work perfectly out of the box, but if you wanted to test it, I'm sure the authors would help you get it to run and value your feedback).

If I want to calculate the gnomAD_AF, which parameter should I add to the nf-core/sarek pipeline?

This is an annotation item (i.e. fixed for each [known] variant, not calculated), and it is included when you have

vep

in your

--tools

. Probably

snpeff

has it as well, but I rarely use that. (All of the VEP results are included in the VCFs in a very-hard-to-read text sausage (the

CQ

INFO

field, so you'll have to extract that, e.g. with the

bcftools +split-vep

plugin.) Note that there are several related paramters, e.g. for each population (Africans, Ashkenazi Jewa, etc.), and depending on the VEP version, this may be split between their WES and WGS data as well (into

gnomADe_AF

and

gnomADg_AF

), so be sure to specify the correct one.