[santiagorevale]

Hi Genomax,

I'm using 10Xs cellranger pipeline.

In this particular configuration, they include both "cell barcode" and "UMI" in Read 1 while the current RNA data is encoded in Read 2. The Index 1 file is used as part of the demultiplexing process because 10X uses 4 Illumina indexes per sample instead of just one.

One additional question: in this scenario, Read 1 is only 26 bp while Read 2 is 75 bp. I want to apply a basic filtering (only maxns) on Read 2, is there a way of doing that? I've tried skipr1=t but that only applies to kmer opperations.

Thanks again.

[GenoMax]

Can you give examples of programs you are referring to? I assume the index read here is not the standard Illumina index but some sort of UMI?

BTW: Generation of index read in a separate file is an atypical application which requires a separate option for the bcl2fastq pipeline.

[santiagorevale]

Hi Brian,

Would it be possible to add to bbduk the capability of filtering also the Illumina Index file? Just to clarify what I mean, Illumina produces the following 3 output files for a typical pair-end single indexed sample:

Read 1 : Sample_S1_L001_R1_001.fastq.gz
Read 2 : Sample_S1_L001_R2_001.fastq.gz
Index 1 : Sample_S1_L001_I1_001.fastq.gz

When I'm filtering using bbduk, I specify in1=read1 in2=read2 out1=filtered.read1 out2=filtered.read2.

What I would like is to be able to add "in3=index1 out3=filtered.index1". I don't think any quality filtering should be applied to this file except only outputting only the sequences that passed the quality filtering for reads 1/2.

Currently, what I'm doing now is extracting the IDs of the filtered file and then filtering the index1 file by these IDs, which is very inefficient.

Many new protocols (mainly for single-cell analysis) are using these 3 Illumina files as input in their pipelines. So I think it could be a good addition.

Thank you very much in advance!

Cheers,
S

[bio_d]

aligning Illumina mate-pair library using BBMAP

Code:

bbmap.sh rcomp=t rcs=f in=read1.fq in2=read2.fq out=mapped.sam

Is this the correct way to align illumina mate pair(5.2Kb)/ long mate pair(10Kb) libraries.

No matter what combination of the flags rcs=t/f rcomp=t/f I use the standard error file shows that "Processing reads in paired-ended mode.". I can't understand how to use the flags correctly because in the UsageGuide it is suggested that one should use requirecorrectstrand=f (rcs=f) and rcomp=t for long mate pair. However, since I am getting insert size of 2611.48 when I am expecting something around 5200bp (predicted 5.2Kb Illumina mate pair library from sequencing center) it is most likely the flags are overridden by the default values which are rcs=t and rcomp=f (which is why I presume bbmap is processing reads in the paired-end mode).

Could you please help.

Best,
D

[Gopo]

Hi Gopo,

Yes, I will add that (as an option). Is that common practice in other variant-callers? Note that callvariants.sh does currently have a "PF" (pass filter) field per sample, but I want to make things as simple and compatible as possible.

Hi Brian,
GATK does output "./." for individuals where a genotype can't be reliably called at a particular SNP (can't remember if freebayes and/or ANGSD does this too). What is the flag for "PF" (pass filter) field per sample?

Is there a flag for filtering --min-alleles 2 --max-alleles 2? I ask because I normally do this with VCFtools, but with big VCF files, this takes a while.

Thank you,
Gopo

[GenoMax]

@Brian: K-mer module in FastQC only tracks 2% of the data (1 in 50 sequences)

[Brian Bushnell]

Hello,

I appear to be getting different results when running the same parameters on the same datasets. I used the below commands:

./bbduk.sh in=R1_001.fastq in2=R2_001.fastq minlen=51 mink=11 ktrim=r tbo tpe K=31 hdist=1 hdist2=1 ref=adapters.fa stats=stats_adaptor.txt out=R1_clean.fq out2=R2_clean.fq outm=ER1_adaptormatch.fq outm2=R2_adaptormatch.fq

In the first instance I had a great improvement in the data quality as reported by fastQC with no problems with kmer content in either read. However, when I had to re-run the same data with the same commands, FastQC now flagged up a warning on my R1 reads on the Kmer content. When I ran it a third time it then flagged up a warning on my R2 reads on the Kmer content.

Each time I run the above commands on the dataset I either get:
-No warning on either set of reads for Kmer content
-A warning on R1 Kmer content but not R2
-A warning on R2 Kmer content but not R1

I see two different kmer warnings, either with one or two overrepresented kmers although what kmers and their position in the read differs each time.

Could anyone shed some light on why this is happening? I've attached screenshots of the two-types of Kmer warning that get flagged by fastqc.

Many thanks for your help as always.

Hi Glfrey,

BBDuk is strictly deterministic on a per-read basis, however it does by default reorder the reads when run multithreaded. You can add the flag "ordered" to keep output reads in the same order as input reads, in which case a dif of the output files between multiple runs with the same input and same parameters should indicate that they are identical.

It is possible that FastQC is subsampling or using the first X reads of the file. In either case, the results could differ slightly if BBDuk is not run with the "ordered" flag (both for FastQC and any process that uses a subset of the reads). Please let me know if using the "ordered" flag does not resolve the issue. The number of reads involved here (as reported by FastQC) seems to be very small. Are the output files of the expected size, or did BBDuk perhaps crash for some reason prior to completion?

[Brian Bushnell]

Hi Brian,
I really like callvariants.sh (simulated di-allelic SNPs and at both low (4x) and moderate (60x) sequencing coverage, filtering to keep only variants with quality score of >=Q27, callvariants.sh is more accurate than using GATK 3.7 with their recommended variant filtering). One thing though that I was wondering about, was whether it would be possible to have an option that sets variants that fail for some individuals (in multisample mode) to a missing genotype call "./." rather than whatever genotype is called by callvariants.sh
Best,
Gopo

Hi Gopo,

Yes, I will add that (as an option). Is that common practice in other variant-callers? Note that callvariants.sh does currently have a "PF" (pass filter) field per sample, but I want to make things as simple and compatible as possible.

[glfrey]

Variable Kmer warning when using bbduk

Hello,

I appear to be getting different results when running the same parameters on the same datasets. I used the below commands:

./bbduk.sh in=R1_001.fastq in2=R2_001.fastq minlen=51 mink=11 ktrim=r tbo tpe K=31 hdist=1 hdist2=1 ref=adapters.fa stats=stats_adaptor.txt out=R1_clean.fq out2=R2_clean.fq outm=ER1_adaptormatch.fq outm2=R2_adaptormatch.fq

In the first instance I had a great improvement in the data quality as reported by fastQC with no problems with kmer content in either read. However, when I had to re-run the same data with the same commands, FastQC now flagged up a warning on my R1 reads on the Kmer content. When I ran it a third time it then flagged up a warning on my R2 reads on the Kmer content.

Each time I run the above commands on the dataset I either get:
-No warning on either set of reads for Kmer content
-A warning on R1 Kmer content but not R2
-A warning on R2 Kmer content but not R1

I see two different kmer warnings, either with one or two overrepresented kmers although what kmers and their position in the read differs each time.

Could anyone shed some light on why this is happening? I've attached screenshots of the two-types of Kmer warning that get flagged by fastqc.

Many thanks for your help as always.

[Gopo]

Hi Brian,
I really like callvariants.sh (simulated di-allelic SNPs and at both low (4x) and moderate (60x) sequencing coverage, filtering to keep only variants with quality score of >=Q27, callvariants.sh is more accurate than using GATK 3.7 with their recommended variant filtering). One thing though that I was wondering about, was whether it would be possible to have an option that sets variants that fail for some individuals (in multisample mode) to a missing genotype call "./." rather than whatever genotype is called by callvariants.sh
Best,
Gopo

[dariober]

I don't like to alter sequence headers by default because that could cause the introduction of non-unique headers

Hi Brain- Thanks for the comments (and for the tools of course!). I see your point but I still think the current default behaviour of bbmap is incorrect when sequence names have spaces because it creates an inconsistency between contig names in the sam header versus the names in the sam records.

Since spaces in names are not allowed by the sam format, I would expect bbmap either to silently trim at the first space (which is not nice...!) or fail straightaway with noting that if you really want to proceed you have to enable the tdr=t option. Just a thought...

[Brian Bushnell]

reformat.sh has an option "underscore" which will change whitespace in sequence headers into underscores, if the extra information is important. Alternatively, as Genomax says, you can use "trimrname". Generally I don't like to alter sequence headers by default because that could cause the introduction of non-unique headers, which would break all tools, rather than just some tools

[GenoMax]

You can choose to truncate the names at the first space during mapping using "trd=t".

Since you already have sam/bam files you can use reformat.sh to take care of the spaces using

trimrname=f             For sam/bam files, trim rname/rnext fields after the first space.

[dariober]

Sorry... another issue.

I have a reference genome where the header lines contain spaces, like this:

>chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38
NNNNN...

After mapping with bbmap I get a sam header where the SN field contains the full fasta header:

cat -vet test.sam | head

@HD^IVN:1.4^ISO:unsorted$
@SQ^ISN:chr1  AC:CM000663.2  gi:568336023  LN:248956422  rl:Chromosome  M5:6aef897c3d6ff0c78aff06ac189178dd  AS:GRCh38^ILN:248956422$

the sam records contain only the first word, i.e. "chr1", "chr2" etc. This causes problems to downstream tools since "chr1 AC:CM000663.2 ..." is not a valid reference name.

I think the printing of the sequence name in the sam header should strip everything after the first blank space?

[dariober]

Hi Dario and Genomax,

I though I had gotten rid of that a while ago, but I guess not - I'll investigate and fix it. It's harmless and due to a race condition for a thread finishing after it was prematurely shut down because enough reads were generated, but it looks scary.

Hi Brain, thanks for feedback. I thought it was harmless the problem is that I have a pipeline (snakemake) that will "refuse" to proceed on receiving this error.