text:
How to normalise Bulk RNA-seq data to account for transcript length, coverage depth, library size and batch effects?
I'm currently working on a project where I'm comparing aggregate measurements (mean, median, etc.) of expression data (RNA-seq) from different groups of genes across various samples with different characteristics (tissue type, health status, etc.). Additionally, the raw counts were collected from several different labs using various techniques.
Since I am conducting between-gene measurements, the data should be normalised to account for differences in transcript length and coverage depth (TPM, RPKM, FPKM). However, I am also interested in comparisons across samples based on tissue type and other factors. Therefore, the data should also be normalised to account for library size (TMM, quantile, etc.), and, as the data were collected from multiple sources, it should be corrected for batch effects.
I have read through many papers but am unsure and confused about how to proceed with the normalisation procedure starting with the raw counts. Can I simply string the methods together, starting with batch effect correction, followed by library size normalisation, and then the within-sample normalisations?
I would appreciate any insights or suggestions on this. Thanks

Please be more specific about your experimental design (e.g. number of samples, factors involved, inputs, desired outputs). General questions don't usually get answered well on Bioinformatics Stack Exchange.

All factors you mention bar the batch effects can be addressed by some sort of normalization (I use a modified version of UQpgQ2). Batch effects can be addressed using ComBat/SVA/ComBat-seq but my team's experience with them has been bad - they introduce artifacts and skew the data quite a bit especially if the batch effect is not pronounced.
Unfortunately, there is no way to do simple comparisons without huge caveats. You're better off using DESeq2 to perform statistically sound comparisons using the batch variable as a covariate.

Have a look through the DESeq2 manual. It goes into a lot of detail about how it works, and how you can use it to do all the things you're asking about:
https://www.bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The standard workflow already has all these things strung together, and works fairly well for the most common differential expression tasks:

Quick start
Here we show the most basic steps for a differential expression analysis. There are a variety of steps upstream of DESeq2 that result in the generation of counts or estimated counts for each sample, which we will discuss in the sections below. This code chunk assumes that you have a count matrix called cts and a table of sample information called coldata. The design indicates how to model the samples, here, that we want to measure the effect of the condition, controlling for batch differences. The two factor variables batch and condition should be columns of coldata.
dds <- DESeqDataSetFromMatrix(countData = cts,
                             colData = coldata,
                             design= ~ batch + condition)
dds <- DESeq(dds)
resultsNames(dds) # lists the coefficients
res <- results(dds, name="condition_trt_vs_untrt")
# or to shrink log fold changes association with condition:
res <- lfcShrink(dds, coef="condition_trt_vs_untrt", type="apeglm")


text:
Discrepancy in Depth of Coverage Estimation between Pileup and samtools depth
I'm encountering a discrepancy in the depth of coverage assessment for specific coordinates when comparing results obtained from pileup files and samtools depth.
I've processed pileup files and extracted coverage information for a given coordinate, resulting in a coverage depth of 8028. However, when using samtools depth on the same coordinate, the coverage depth is reported as 88565.
My analysis involves the following steps:
Alignment of paired-end reads using bwa mem.
Generation of pileup files and extraction of coverage information.
Utilizing samtools depth to calculate coverage depth for the same coordinates.
Visual inspection in IGV for validation, where I observe that the number of reads aligning to the region matches the count obtained from samtools depth.
I've noticed that in this region, there are no structural variations. The versions of samtools I'm using are samtools 1.19.2.
Could you please help me understand the potential reasons for such big discrepancies? What parameters should I review, and could the differences be attributed to any specific factors?
Thank you for your assistance!

Please clarify your specific problem or provide additional details to highlight exactly what you need. As it's currently written, it's hard to tell exactly what you're asking.

How are you generating the pileup files? Can you please show the command line steps you use in your analysis to generate these pileup files & depth calculations?

This is purely a guess based on what you write given that you are not showing the exact command for samtools etc.
I suspect that you are allowing sequences with MAPQ==0 (supplementary alignments) into your depth calculation. A brief google search tells me that this issue has been previously observed with samtools depth.
In principle, you should be able to see this directly in IGV, by controlling how supplementary alignments are displayed.
This is something that you can directly control with the -Q flag, to restrict it to uniquely mapping reads (where you set the flag to e.g. -Q 1 or higher).
I'd suggest trying to play with samtools depth parameters and seeing if you can bring it more into alignment.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/

sort: votes


text:
How to classify TF motifs by family
I have a FIMO output, the input sequences were putative promoter sequences of several species. I want to graph the positions of these motifs along a horizontal graph, but I want to filter the data by clustering the motifs by family, e.g. all the zinc finger motifs being represented by one color, all the helix loop helix motifs being represented by one color, etc. How would I go about clustering the raw .tsv file?

Please explain what FIMO is, and show the first few lines from your raw .tsv file. What language / tool would you like to use for clustering? Is this something you'd do directly in the FIMO tool, or externally using another program?

JohnDoe23, are you looking for an output like the one shown here: https://bioinformatics.stackexchange.com/a/20132/3967

Comment: I don't understand what you mean by "I want to graph the positions of these motifs along a horizontal graph"? Are you talking about a Manhatten plot? Or something else? And I'm not sure about "filter the data by clustering the motifs by family". What type of clustering do you want to perform?
Without further information/description it's difficult to understand what you're trying to do. If possible, please edit your question with an example of what your desired output. Also please indicate the language you're using for your visualisation, e.g. R, python, other.

Answer: I may have completely misunderstood, but perhaps this approach could be adapted to suit. Using the example FIMO output file from https://meme-suite.org/meme/doc/examples/fimo_example_output_files/fimo.html?man_type=web:
library(tidyverse)
library(qqman)
#>
#> For example usage please run: vignette('qqman')
#>
#> Citation appreciated but not required:
#> Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731
#>

fimo_df <- read_tsv("~/Desktop/fimo.tsv", comment = "#") %>%
  rename("SNP" = motif_id, "CHR" = sequence_name, "BP" = start,
         "P" = `p-value`, "zscore" = `q-value`) %>%
  mutate(CHR = case_when(str_detect(CHR, "X") ~ "23",
                         str_detect(CHR, "Y") ~ "24",
                         TRUE ~ CHR)) %>%
  mutate(CHR = parse_number(CHR)) %>%
  select(-c(stop, strand, matched_sequence))
#> Rows: 1494 Columns: 10
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: "\t"
#> chr (5): motif_id, motif_alt_id, sequence_name, strand, matched_sequence
#> dbl (5): start, stop, score, p-value, q-value
#>
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

manhattan(fimo_df, p = "zscore", logp = TRUE, ylab = "q-score",
          genomewideline = FALSE, suggestiveline = FALSE,
          main = "Manhattan plot of TF motifs",
          chrlabs = c(1:19, "X", "Y"))

If you want to colour the points by TF motif (or by 'cluster'), perhaps ggplot would be appropriate, e.g.
library(gghighlight)
axis_set <- fimo_df |>
  group_by(CHR) |>
  summarize(center = mean(BP, na.rm = TRUE))

ylim <- fimo_df |>
  filter(P == min(P)) |>
  mutate(ylim = abs(floor(log10(P))) + 2) |>
  pull(ylim)

ggplot(fimo_df, aes(x = BP, y = -log10(P),
                    color = as.factor(CHR),
                    size = -log10(P))) +
  geom_point(alpha = 0.75) +
  # scale_x_continuous(label = sort(axis_set$CHR), breaks = sort(axis_set$center)) +
  scale_y_continuous(expand = c(0, 0), limits = c(3.95, ylim)) +
  scale_color_manual(values = rep(
    c("#276FBF", "#183059"),
    unique(length(axis_set$CHR))
  )) +
  scale_size_continuous(range = c(0.5, 3)) +
  theme_minimal() +
  theme(
    legend.position = "none",
    panel.grid.major.x = element_blank(),
    panel.grid.minor.x = element_blank(),
  ) +
  gghighlight(motif_alt_id == "TBP")
#> Warning: Tried to calculate with group_by(), but the calculation failed.
#> Falling back to ungrouped filter operation...

Created on 2024-02-05 with reprex v2.1.0
Does that help at all?

How do you know the "family" to which your motifs belong?  If the source for the motifs that FIMO matched is, say, JASPAR Vertebrata CORE, then you might reasonably decide to use JASPAR's "structural class" as provided by their pre-computed "Annotation table" in JASPAR CORE Vertebrates clustering.
Assuming you are using JASPAR and you do make this decision...
Your "raw .tsv file" is in FIMO Results TSV Format.  "Joining" it with a downloaded copy of the "Annotation table" serves to annotate it with additional columns, having the effect you seek of clustering the FIMO results on the "structural class" of each matched motif.  Joining can be implemented following Running queries directly against CSV or JSON using the sqlite-utils command-line utility.
How you would then plot these features and encode the class as a color requires more detailed description what the plot is supposed to represent or look like.  For instance, do you "want to graph the positions of these motifs along a horizontal graph" with coordinates taken with respect to the TSS of the transcript from which the promoter is taken?

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
WGCNA: Extract genes within the different modules to create a dataframe with two columns
My inquiry pertains to WGCNA. After identifying the modules, I would like to extract genes within the different modules to create a dataframe with two columns: Gene_a and Gene_b. The values in these columns should represent co-expressed genes. An illustrative example is as follows:

In this example, the gene XX62730x is co-expressed with XCGR1Y. And YYYYY is co-expressed with ZZZZZZ.
However, I am facing difficulties as the official WGCNA page is no longer accessible.
I would greatly appreciate any help you can provide. Thank you in advance!
Best,
Nana

The archive.org copy might help while you're trying to figure out the answer: https://web.archive.org/web/20231002045159/https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA/

@RamRS thanks for the link!

Unfortunately, the links on that link  do not work.

I tested a couple and they worked - you need to learn to navigate archive.org to be able to leverage those links the best way.

U can follow WGCNA R packages

You can get this information, more or less, using the functions exportNetworkToVisANT and exportNetworkToCytoscape. You will find some extra columns that you can ignore or delete from the resulting table. We are still working with UCLA to restore the WGCNA page including the tutorials, but the tutorials can temporarily be accessed using this dropbox link.
To use either one of the two functions, you will need the network matrix for which we usually use TOM. You will probably want to set a threshold for connections to be included, otherwise you will get a fully connected network and the resulting table will (usually) be huge.
Hope this helps!

thank you very much, that really helped me


text:
How to build hisatgenotype extra index files (.link, .locus, etc)
I'm trying to run hisatgenotype on a custom genome, originally specified in fasta format. I'm able to use hisat2-build to build the .1.ht2, .....8.ht2 files, however I've not been able to figure out how to get the other files that hisatgenotype requires. I've looked into hisatgenotype_toolkit build-genome but it errors out due to problems with argument parsing. Is this even possible, or do I have to use the genomes/indices downloaded by the setup.sh script?
edit:
The error I get is the same as this issue: https://github.com/DaehwanKimLab/hisat-genotype/issues/69
% hisatgenotype_toolkit build-genome --base hla
Traceback (most recent call last):
  File "/Users/ted/miniforge-x86/envs/hisat/hisat-genotype/hisatgenotype_tools/hisatgenotype_build_genome.py", line 546, in <module>
    typing_common.clone_hisatgenotype_database()
  File "/Users/ted/miniforge-x86/envs/hisat/hisat-genotype/hisatgenotype_modules/hisatgenotype_typing_common.py", line 52, in wrapper
    func(*args, **kwargs)
TypeError: clone_hisatgenotype_database() missing 1 required positional argument: 'destination'
Script exited with error 1

If you try to run with an extra positional argument it fails too
% hisatgenotype_toolkit build-genome --base hla indices
usage: hisatgenotype_build_genome.py [-h] [--base BASE_FNAME]
                                     [--locus-list LOCUS_LIST]
                                     [--inter-gap INTER_GAP]
                                     [--intra-gap INTRA_GAP]
                                     [--aligner ALIGNER] [--linear-index]
                                     [--commonvar] [--clinvar] [-p THREADS]
                                     [-v]
hisatgenotype_build_genome.py: error: unrecognized arguments: indices
Script exited with error 2

And similar problems if you add --destination indices

"but it errors out due to problems with argument passing" - show us your exact command and the error being displayed.


text:
limma linear model with interaction
I want to run some epigenome wide association studies based on city and pollutant type, could you please suggest if I can get DNA methylation sites (CpG) association based on pollutant type for each city. like below I am trying to find Pb toxic elements association site in four different cities which has varying level of Pb.
#model matrix
var<-model.matrix(~log2Pb*City + sex + Age +CD8T +CD4T +NK + Bcell +Mono, data=targets2)

fit<-lmFit(mval,var)

fit2<-eBayes(fit,trend=TRUE, robust=TRUE)

probe<-topTable(fit2,adjust="BH",coef=2,num=Inf)

try to run interaction model using limma but not sure if this is the right approach


text:
In medical and public health care, should we make use of "clinical" signficance testing (non-parametric tests) or parametric tests of significance?"
In certain cases, traditional significance test may not be able to identify that effect of a treatment or a drug is significant. However,clinical significance test may suggest a significant effect. Does it mean that this test is a better way of interpreting  effect-size or correlation coefficient when we are facing imperfect measurements or data. Imperfect measurements could produce a wrong conclusion based on traditional statistical test?

Your question is not clear. Can you give us an example where the significance test indicates a lack of significance (with good experimental design) but other evidence says otherwise?

I'm sorry but this is just too broad a question. If you [edit] and give a specific example we might be able to help if it is about bioinformatics, but as it stands, we cannot help you assess the results of a statistical test you haven't shown.

@RamRS please have a look. https://bioinformatics.stackexchange.com/questions/11043/difference-between-clinical-significance-and-statistical-significance

Thank you for the pointer, @SubhashC.Davar

@Ram RS Meta Analysis: please read my paper - An Improved Version of Hunter, Schmidt and Jackson Framework.


text:
Loop and TAD calling from HiC reads using HOMER
I have HiC reads from an HiC experiment and I am using HOMER to call loops and TADs. The numbers I am getting are like:
Total TADs: 5024; Total Loops : 254
As multiple loops are expected in a TAD so numbers of loops should be higher than the number of TADs. That's why I am not sure how to justify the numbers I am getting.
I made sure I am not doing any mistake in the HOMER pipeline of loop/TAD calling and I followed all the steps in the tutorial here: http://homer.ucsd.edu/homer/interactions2/index.html Unfortunately I don't have any other loop calling program installed right now and HOMER is easy to install and run.
If someone working with HiC data analysis and TAD/loop calling from HiC reads can please comment on this?
Thanks
Following is the set of commands I used (based on HOMER tutorial here http://homer.ucsd.edu/homer/interactions2/index.html) to get loop/TAD calls from HiC reads.
#Trimming/processing of HiC fastq files:
homerTools trim -3 GATC -mis 0 -matchStart 20 -min 20 Hu_concat_1.fastq
homerTools trim -3 GATC -mis 0 -matchStart 20 -min 20 Hu_concat_2.fastq
#Aligning the two paired-end fastq files to get sam:
bowtie2 -p 20 -x hg38.fa -U Hu_concat_1.fastq.trimmed > Hu_concat_1.sam
bowtie2 -p 20 -x hg38.fa -U Hu_concat_2.fastq.trimmed > Hu_concat_2.sam
#Making tag directory (the HOMER makeTagDirectory program is used to process the alignment file into a HOMER-style tag directory which is used for subsequent analysis steps)
makeTagDirectory Hu_Hic_TagDir/ Hu_concat_1.sam,Hu_concat_2.sam -tbp 1
#Calling loops and TADs
findTADsAndLoops.pl find Hu_Hic_TagDir/ -cpu 10 -res 10000 -window 15000 -genome hg38.fa
The whole process ran smoothly without errors. It is just that in the output, the number of loops I am getting are far fewer than no. of TADs (Total TADs: 5024; Total Loops : 254) which I think is not biologically correct. If someone using HOMER for calling loops and TADs from HiC reads can justify this output, that would be really helpful.


text:
Fastest way to convert GTF or GFF file to a BED file?
Given a GTF/GFF3 annotation file, what is the fastest way to convert it into a BED file?

Hi @AlejandroGonzalesIrribarren. Thanks for your response. For future reference please provide some context to this question, e.g. the genome (humans), computer resources and what have your tried? Again thanks for posting back. Conda is cool for getting stuff working quickly, but not necessarily performant .... Easybuild is much better in this particularly question, i.e. fastest.

Yes, please at least provide example input and expected output so the question can be clear.

gxf2bed
The fastest G{F,T}F-to-BED converter around the block!
translates
chr27 gxf2bed gene 17266470 17285418 . + . gene_id "ENSG00000151743";
chr27 gxf2bed transcript 17266470 17281218 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8";
chr27 gxf2bed exon 17266470 17266572 . + . gene_id "ENSG00000151743"; transcript_id "ENST00000541931.8"; exon_number "1"; exon_id "ENST00000541931.8.1";

into
chr27 17266469 17281218 ENST00000541931.8 1000 + 17266469 17281218 0,0,200 2 103,74, 0,14675,

before your eyes blink!
Converts

Homo sapiens GRCh38 GENCODE 44 (252,835 transcripts) in 3.34 seconds.
Mus musculus GRCm39 GENCODE 44 (149,547 transcritps) in 2.11 seconds.
Canis lupus familiaris ROS_Cfam_1.0 Ensembl 110 (55,335 transcripts) in 0.99 seconds.
Gallus galus bGalGal1 Ensembl 110 (72,689 transcripts) in 1.14 seconds.

Usage
Usage: gxf2bed[EXE] --input/-i <GTF/GFF> --output/-o <BED>

Arguments:
    --input/-i <GTF/GFF>: a .gtf/.gff file
    --output/-o <BED>: path to output .bed file

Options:
    --help: print help
    --version: print version
    --threads/-t: number of threads (default: max ncpus)

Installation
to install gxf2bed on your system follow this steps:

get rust: curl https://sh.rustup.rs -sSf | sh on unix, or go here for other options
run cargo install gxf2bed (make sure ~/.cargo/bin is in your $PATH before running it)
use gxf2bed with the required arguments
enjoy!

Build
to build gxf2bed from this repo, do:

get rust (as described above)
run git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed
run cargo run --release -- -i <GTF/GFF> -o <BED>

Container image
to build the development container image:

run git clone https://github.com/alejandrogzi/gxf2bed.git && cd gxf2bed
initialize docker with start docker or systemctl start docker
build the image docker image build --tag gxf2bed .
run docker run --rm -v "[dir_where_your_gtf_is]:/dir" gxf2bed -i /dir/<GTF/GFF> -o /dir/<BED>

Conda
to use gxf2bed through Conda just:

conda install gxf2bed -c bioconda or conda create -n gxf2bed -c bioconda gxf2bed

Benchmark + FAQ

If you google "gtf to bed" or "gff to bed" you'll find some posts about recommended tools. Part of these tools are already deprecated or too difficult to run (even in Linux), due to poor maintenance. This project was conceived to provide an easy way to convert GTF/GFF files to BED structures and finish a set of high-performance converters between gene model formats in Rust (bed2gtf, bed2gff and now gxf2bed).
To test the efficiency of gxf2bed, I took 4 GTF/GFF-to-BED converters and ran them with the current GRCh38 GENCODE annotation (v.44), that has ~250,000 transcripts [the values showed here are the mean of 5 consecutive runs]. Briefly, the competitors are:

UCSC's utils: gtfToGenePred | genePredToBed & gffToGenePred | genePredToBed (taken from here)
Signal & Brown: gtf2bed.py & gff32gtf.py && gtf2bed.py (taken from here)
ea-utils: gtf2bed.pl & gff2gtf.pl | gtf2bed.pl (taken from here)
BEDOPS: gtf2bed &  gtf2bed (taken from here)

Besides the easy way to make gxf2bed run and the GFF-GTF channel that allows centralize both in a single converter, this tool showed a significant difference in computation time against the other tools in each one of the two formats. Even using a combination of different number of threads, times were practically the same 3.5s +/- 0.4s and the differences were maintained. Is important to note that GFF3 files took a lot more time than its GTF counterpart. This could be due to the need of chaining two different programs 1) GFF-to-GTF and 2) GTF-to-BED.
Taken together, gxf2bed offers an easier and faster way to convert GTF/GFF files to BED files. This tool could save you at least x2-3 times the time you used to spend using other tools to convert GTF files and x5 times if you want to convert GFF3 files!

Disclaimer: I am the author of the gxf2bed tool.

How do you deal with regular GTFs that have a variable number of key-value pairs in the ninth column?

It does not matter how many number of key-value pairs are in the ninth column. The only thing this tool assumes is that at least "gene_id" and "transcript_id" should be present in each line. This is nothing crazy to ask since is very common practice.

So what options did you use to make the conversion you show here? Why were the gene and exon lines ignored and only the transcript kept?

Exactly what @terdon says - there doesn't seem to be any way to pick the features to extract. Is your tool faster than a simple awk in that case? And are you sure you're comparing apples to apples when your tool makes these assumptions?

Exon lines are not ignored. The tool basically parses a gtf and extracts only the needed features to build a BED12 (which is basically a transcript-based annotation file). The results of gxf2bed are exactly the same obtained by UCSC's binaries and the other tools as well. gxf2bed even goes a step further sorting the BED12 by chr and pos, which does not seem to be accurate for some of the tools evaluated.
If you are confident that an awk command could build a correct BED12, wonderful! This tool provides another opportunity to perform this conversion in less time that the available tools now.

If you are skeptical about how the results of gxf2bed, you can try it, though. Just run the command as explained in this post and compare the results yourself. To elaborate a little bit more on the last lines of the past comment: re-writting this functionality in Rust now allows 1) a more straightforward installation, 2) faster times, 3) easy pipeline implementation, 4) better performance. All this is open source, so feel free to review the code if you still do not feel that comfortable! :)

"a BED12 (which is basically a transcript-based annotation file)" Please show your source for this information. BED files are lists of contigs+positions with optional annotations, they could be genes, transcripts or even random nucleotides. You are assuming the transcript part.

gtf2bed, for example, accepts an attribute_key and is thus flexible. gtf2genePred does assume gene as the feature and it does provide some flexibility as to the identifier within that constraint, but it calls itself genePred, not featurePred. Sure, gtf2bed.pl seems to default to transcript_id but that's the only tool you can compare to. I appreciate your rust-based contributions to the field but please be more mindful with your claims.


text:
Finding homologous regions in multiple whole genomes
Right now I have six genomes that I want to compare and identify homologous regions in the genomes. I have run nucmer, show-coords and obtained the output files. An example is shown below with Genome 1 vs Genome 2. More files go through Genome 1 vs Genome 3 until all are compared.
I was trying to use bedtools merge after concatenating the multiple output files from show-coords command to find the homologous regions. However, I am not quite sure on the command to use. Any help is appreciated.
Edit: added commands. I just want to know what regions are homologous across the six genomes (99/100% IDY). I don't know what I am supposed to merge to find this out. I was reading and I believe bedtools should be but I don't know. I am usually a wetlab but I try to do a bit of bioinformatics.
 # Loop through all pairs of genomes
for i in {1..6}; do
    for j in $(seq $((i + 1)) 6); do
        genome1="genome${i}.fasta"
        genome2="genome${j}.fasta"

        output_prefix="output_${i}_${j}"

        # Run NUCmer
        nucmer -p "${output_prefix}" "$genome1" "$genome2"

        # Show coordinates
        show-coords -r -l -c -T "${output_prefix}.delta" > "${output_prefix}_coords.txt"
    done
done

After using the above I then concatenate the files with
cat output_*_coords.txt > combined_coords.txt

This gives me something that looks like this below.
[S1]    [E1]    [S2]    [E2]    [LEN 1] [LEN 2] [% IDY] [LEN R] [LEN Q] [COV R] [COV Q] [TAGS]
1       14347   1       14347   14347   14347   99.02   5064019 5004926 0.28    0.29    Genome1      Genome2
14474   32175   14350   32051   17702   17702   99.33   5064019 5004926 0.35    0.35    Genome1      Genome2
16541   22042   4698474 4692975 5502    5500    99.64   5064019 5004926 0.11    0.11    Genome1      Genome2
16548   22042   3728461 3722968 5495    5494    99.84   5064019 5004926 0.11    0.11    Genome1      Genome2
16583   22041   3904238 3898781 5459    5458    99.65   5064019 5004926 0.11    0.11    Genome1      Genome2
18727   22097   4426876 4423506 3371    3371    99.73   5064019 5004926 0.07    0.07    Genome1      Genome2

Please [edit] your question and show us the exact commands you used. All of them. Including how you concatenated the files. We need to know exactly what each of the columns is and then we need to see or at least understand your expected output. You only say you want to use bedtools merge, but you don't tell us what you are merging or what you want to do. When should a line be merged? When it is 100% identical? More than N% identical? Any size limits? How much identity across genomes? Please explain what you need in more detail.

According to show cords the format you are getting contains columns which you can use to obtain BED format with relevant information.
You need following columns:
TAGS_genome1
[S1]    Start of the alignment region in the reference sequence.
[E1]    End of the alignment region in the reference sequence.
TAGS_genome2
[% IDY]

I would start with not merging all coordinates but getting bed files for each pairwise genome_a__vs__genome_b.bed etc. and start from there, taking a subset of 3 genomes first (say a, b,c).
You can find the regions shared by all three bedtools intersect:
bedtools intersect -a a_vs_b.bed -b a_vs_c.bed

The issue is that you will have many to many relationships, since region_1 from b can have several similar regions in b, and each one of these (== regions in b) also can be similar to a number of regions in a. Depending on what you are looking for and how dis-similar the genomes are you may end up with reporting repetitive elements or if you filter hard missing i.e. exons.

Thanks @darked89. That's exactly what I needed, those 5 columns. I didn't consider that there are too many relationships.
I will then continue to do comparisons. I think will give me the answer that I want.
bedtools intersect -a a_vs_b.bed -b a_vs_c.bed > a_vs_b_vs_c.bed
bedtools intersect -a a_vs_b_vs_c.bed vs a_vs_d_vs_e.bed > all_compared.bed


text:
split fastq file containing a sequence block at different locations
I have some fastq files (obtained from nanopore sequencing) that contain reads that can be of either of these 5 forms:

a known CDS with 3'UTR:

CDS----------------------- (seq1: original sequence)

the same cds with a block of 50bp represented by ||||| (that substitutes the 3'UTR at 4 different locations)

CDS|||||------------------ (sub1: substitution 1)
CDS-----|||||------------- (sub2: substitution 2)
CDS----------|||||-------- (sub3: substitution 3)
CDS---------------|||||--- (sub4: substitution 4)

I need to split my fastq files in 5 smaller files clustering the different types of sequences (seq1, sub1, sub2, sub3 and sub4). The problem is that those sequences all share some segments and I don't know which method I should apply. I tried finding the best match by building a BLAST database with the 5 sequences but as it is local alignments, I get HSPs for each sequences and cannot discriminate between them.
Does anyone has an idea how I could achieve that?
Any help would be amazing.
edit
Here are 5 sequences as I expect them (i.e., fasta reference)
>seq1 (this one does not have the 50bp block)
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub1
_GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub2
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub3
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACT _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ TTGTCTCTTTTGGTTTATTTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA
>sub4
ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTA _GTAAGTGAGGGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACA_ GATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAA

examples (50bp block if present is in lowercase here)

this read should be associated with sub3 because of the flanking sequences (...CTACAACT in 5' and TTGTCTCTTTT... in 3' of the 50 bp block)

@e213fdd5-bbfb-4dfe-9b03-29eea256da0d
TGTGTACTTCGTTCAGTTACGTATTACTAGTTATTGAGTGTCTTTGTGTTTCTGTTGGTGCGTCTTCGCACAAGGCTAATCTTACTCTCCTCTCCAATGTACTGCTTCTATGTTCATCTCAGCCGAATTCAATAAGGAGAAGAACTTTCACTAGGTTTCCCCATTCTGTTCGATTAATTGGTGATGTTAGTAAATTTCAATTTTCTGTCGGTGGAAGGTGAAGGTGATTGCTGCTGGCTTTTCATTGTTTACTTTTGGCTTCTGTTTTGTAGTAACGCATCACTACTTTCTCTTATAGTGTTCAATGCTTTCCAAGTCTGGAATCATGAAGCAGCGGCTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCAGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACGTGCTGAAGTCAAGTTTGGGGAGACACCCTCGTCAATGAGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACTACAACTCCCACAACGTATACATCATAGCCGACAAAACCCAAAAGACGGCATCAAGAAACCAACTTCAGACCTTACCATCAGGCGGCGGCCATCACACAGTCATTGTCAACAAATGTAATTGGCGGTGGCTATTCTTTTTACCCAGACAACCATTACCTTGTCCACACAGTCACACTGTGAAAAGATCCCAACTTTGAAAGAGAGGCCACTGGTCTTCTTGAGTTTGTAACAACTGCTGGGATTACACATGGCATGGATGAACTATACAAACATGACAGACTCTAAAGCCTGCAGTTCCTGAAGAGATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAGTGACGATGGAACCTGTGTCTACAACTgtaagtgagggagtctcggactcaattagaggctcttctttcacaTTGTCTCTTTTGGTTTATTTTCTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTTTTGTTTTCCTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATGCAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAAACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACTG
+
$&&&%234441124--585545;>B99968650-++*()'&'(;D540488?=888,++*+*)*$$$##&%$$$%)*****('%$$'1/+++-,('))'%$$%$$$&)&%%%&(*'&),/'''&&'06872..*&%)),,3/49:;<?A773)('()/066+'&&)('(+&%$&''$$&&()*--888)('()'')')'''++*+)&')&$$$%()'')*''&&))&'&$&(&$#$%&%&&&%%$%&..&%%%%&&%$%&'&$$'-)(%%%$&(*,+****&'''('+17<{..+*('((*+)*&&&'**.*%$&%&&%$$$%&&'&&(+'(&&&%%&&($$%&&'+55536988;;74632+(*)''%6<=(((&*''(68998===?ABBBB=B?<;;96667>>?>>?@?A=-+.7556<@>A<0//002(()-989::::;<<>>;:5++++->;::;<?F:7<()988;;<?>9;9)36:5-963/002B@A@@?554543,,9<=<@?@>777461166.../B88<;<@=<:888;;>42(''%%(%%+*('&%%&&**)(()*+020*'('&$$$'('''&%'((%(&$##$&)*()&%%%$$#$%'&(&'%$$%%&'%%%$%%%$&('$##$&&&&$$$%$%$$%&%%&040(((&'.023266,,,++0)++-.-*,,0+&&'$&''''%%'&-+,0((+))*+*(%&)+00.3)''%'&'''(*(%%)++++,--+.....0'(()7:<=;=<<<<<0.<;9:;<@;;;:<<=;2+)++,1212611)())*:1++*'&&%((&&056/,454))*''''(3346;<;;;;;000.0/043546??=;<*)+64335663455665333788:78779<>9755446::<><198662*+7899722248867>?8++4))+6779<<BA><9;<;:;;{<?BB-,,-13444/154///933{3==@B:9989::::::;;<<=:3335=EBB;:*)..'-:69?@CKKKDC<<?>:)()((012.,,+,55;=?>@>:98:78544344*)%$$$&)))))))))(((((''''&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&')+.1866677978:;;:::;;;99::988889653456:::9889;::8899.,+0668:;:::8436735//:632/-,+**+013459;9332/+($

this read should be associated with sub4 because of the flanking sequences (...TTTTACTCCTA in 5' even if a C is missing, and GATCTTGTCCT... in 3')

@107dfd6d-2101-4fd8-aeb5-2c432e176d72
ATGTACTTCGTTCAGTTACGTATTGCTATCGCCTACCGTGACAAGAAAGTTGTCGGTGTCTTTGTGTTTCTGTTGGTGCTGATATTGCATGAAGACTAATCTTTTCTCTTTCTCATCTTTTCACTTCTCCTATCATTATCCTCGACCGAATTCAGTAAAGGAGAAGAACTTTTCACTGGGTTGTCCCAATTCTTGTTGAATTAGATGGTGATGTTAATGGGTACAAATTTTCTGTCAGTGGAAGAGGGTGAAGGTGATGCAACATACGGAAAACTTACCCTTAAATTTATTTGCACTACTGGAAAACTACCTGTTCCATGGCCAACACTTGTCACTACTTTCTCTTATGGTGTTCAATGCTTTCAAGATACCAGATCATATGAAGCGGCACGACTTCTTCAAGAGCGCCATGCCTGAGGGATACGTGCGGGAGAGGACCATCTTCTTCAAGGACGACGGGAACTACAAGACACGTGCTGAAGTCAAGTTTGAGGAGACACCCTCGTCAACAGGATCGAGCTTAAGGGAATCGATTTCAAGGAGGACGGAAACATCCTCGGCCACAAGTTGGAATACAACCGCAACTCCCACAACGTATACATCATGGCCGACAAGCAAAAGAACGGCATCAAAGCCAACTTCAAGACCCGCCACAACATCGAAGACGGCGGCGTGCAACTCGCTGATCATTATCAACAAAATACTCCAATTGGCGATGGCCTGTCCTTTTACCAGACAACCATTACCTGTCCACACAATCTGCCCTTTCGAAAGATCCCAACGAAAAATTTTTGCATGGTCCTTCTTGAGTTTGTAACGGCTGCTGGGATTACATGGCATGGATGAAGAACTGCAACACAAACACTTGACGAACTCTAAACCTGCCAGTTCCTGAAGAAATTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGCTAGAAAGTGACGATGGAACCTGTGTCACCAACTTTTGAGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCTAgtaagtgagggagtctcggactcaattagaggctcttctttcacaGATCTTGTCCTCTTTTGTTCCCTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTCATAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAACTGTAGGCACCATCAATGAAGATAGAGCGACAGGCAAGTCACAAAGACACCGACAACTTTCTTGTCACGGTAGGCGATAGCAATACGTAACT
+
-0,+))),:::<??===>=<<<?@>=<<=?;:889:<;<<=;;;;CD@?:9:A>9999=999DBAA@<BE=3220111<B@@AA@A?>>>=>>;:9999:9;<E9A<;?????C8788;==@?;;;==9989@A6677?=7788+**+368>B87::==@<99;<<<<<<??>@=:9/.))*698.*,/8A=<==>>=?BAB=::::;;;<?@>@B@@?:888:>?@CCCEDB>??>===::.+)()2157=7:99<B?>???@??=;::=DGC=<6?;84568>AB@CCBCBAA====AAACEI{AA=<::;<?@=<<>=>@AAA@>>>>?<?====A?==<=@>??@>?==??CA@AAA@AB=757:965(*.;<;<::?AAA@>=>::::<=?@BCC5336899=?>>?@A??42220;9:89.-***+77;97<;:;<<??BC{878<:9:<:2112<<>===@=<;;<>???877/,.-.-///0,++,)299:=;:;;=@@<;=7778?=::;;<;=<<C>AB<BBBCEDCDC?A>;;9:655:DAA@?>>>=?=><;<<@=>><=B=@=:870++,,;=;>><<==CBA@?>=>>=667<;;7778B>B?B<<==<99:::4:659;6642212:.-<=??=6556;??@@>::;:==;8887899==;2143235<<<?>==?CEBD?>>>ADB=720/.02334<00//06568>:888:711*)-776'''&06:<=><><866;{C@A?>8897=>EK558568552-179664)(()),,--,()))3:;:<BAAA>=>>>><<>>*))+/245999;>>=96;?B9888<<77-,,,,&&%%%&&'+6;A;2111--,,.:;><<ACCA@@?=9;;;>=>==@930../+)(())07>>8<><;:;9238<>=>998878:7700016545**1236434689<9<<;.--0-,,,&'():78862237@???><;;;;<<=@>@??@@===;<=>>;:8;)(67511223/.---.2001023:;21<?=**6223;;>?AACGA=73366:;<=:766;ACC=<;86000?8=@DHEC@@=;;<9999<<==9><==8889>A>=;30266777:972+***)))'&%%%%%$$$$$$######"""""""""""""""""""""""""""""""""""""""""#&')+-1122/...46688<<=0('%%&258899;92136:999;:::9:;;--2-.343668=<?@EF????==:7668;6667;><<;;>;;8851.*

@terdon, I added two examples. I cannot provide all cases, as I actually have 16 of them and i'm having a hard time doing it by hand. I still hope that it is clearer now.

You can try seqkit locate with say start or end of CDS sequence. This will give you a tabulated output with locations of your CDS query in each read.
Since you expect that the CDS sequence should be located at one of the five different locations inside your reads you can figure out to which one of the five a given read belongs.
You may also use clumpify.sh from BBMap to cluster reads prior to searching for the CDS locations. The idea is that maybe the five groups will cluster in fairly well separated blocks, but this is fairly speculative.
edit
Since the alignment of your sequences looks like this:
CLUSTAL 2.1 multiple sequence alignment

sub1            NNNNNNNNNNGTAAGTGAGGGAGTCTCGGACTCAAT---TAGAGGCTCTTCTTTCACANN
sub2            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTNNNNNNNNNNGTA
seq1            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
sub3            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGGC
sub4            --ACCTGCCAGTTCCTGAAGGGTGTACTGATCCTGTGGCTGAAAACTTTGATCCAACGG-
                          **   *** **     * **  *  *   *  *  **

sub1            NNNNNNNNTTGATCCAACGGCTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub2            AGTGAGGGAGTCTCGGACT-CAATTAGAGGCTCTTCTTTCAC---ANNNNNNNNNNTTTG
seq1            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
sub3            TAGAAGTGACGATGGAACCTGTGTCTACAACTNNNNNNNNNNGTAAGTGAGGGAGTCTCG
sub4            --------------------CTAGAAGTGACGATGGAACCTG---TGTCTACAACTTTTG
                                              *                          * *

sub1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub2            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
seq1            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTATTGTCTCTTT-----TGGT
sub3            GACTCAATTAGA--GGCTCTTCTTTCACANNNNNNNNNN--TTGTCTCTTT-----TGGT
sub4            AGCAATATTATCCTGCTTATTAATTTGCTGTTTTACTCCTANNNNNNNNNNGTAAGTGAG
                  *   ****    *  * **  **  *                            **

sub1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub2            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
seq1            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub3            TTATTTTT---------------CTCCTTTGTGTAATTGTGGATTGGATCTTGTCCTCTT
sub4            GGAGTCTCGGACTCAATTAGAGGCTCTTCTTTCACANNNNNNNNNNGATCTTGTCCTCTT
                  * * *                *** * * *   *          **************

sub1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub2            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
seq1            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub3            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
sub4            TTGTTCCCTTTTTTTTTTTTTATGATGTACAACACATTGGTAATTTAAAATTGCCTTGTC
                ************************************************************

sub1            ATAAA
sub2            ATAAA
seq1            ATAAA
sub3            ATAAA
sub4            ATAAA
                *****

instead of CDS end/start you have to select regions around your deletions with say 15-20bp flanks on each side. Put these in fasta and run locate.
caveat
running with zero mismatches will miss some reads.
edit 2
The command:
zcat input.fq.gz | seqkit locate --threads 8 \
--max-mismatch 0 \
--pattern-file fasta_with_5_frags_extracted_from_algnt.fa > result.tsv

Some trial and error is required to select good fragments matching exactly one of the five groups plus increasing max-mismatch if due to sequencing errors a lot of reads were dropped.
To speed up IO one can compress result.tsv on the fly, i.e. using pigz or zstd.
In the next step we need to create five separate files listing read names containing a given fragment:
grep frag_name_1 result.tsv | cut -f 1 > frag_name_1.read_names.txt

seqkit grep --pattern-file frag_name_1.read_names.txt input.fq.gz -o frag_name_1.reads.fq.gz

# repeat five times with different pattern file and changing the output

Thanks for all the tips. I actually think that this method can work.
I have tried to use seqkit locate to find the 50bp block and ran it a second time to locate the end of the CDS. Then, if I compute the distance between the two, I get distances roughly multiple of 50, which allows me to distinguish the sequences group I am after. I am going to try and make something quite clean and then I will post here my solution!


text:
How to calculate distance between two regions from bed like dataframe in R?
I'm wondering how can I calculate the distance between two regions from the bed (chr, start, end) alike dataframe in R. Precisely, I want to subtract end of r from start of r+1 (r=row). Also, the calculation will repeat similarly for each chromosome.
head(G3_dat[1:3])

# A tibble: 6 × 3
  chrom   start     end
  <chr>   <int>   <int>
1 chr1  1330426 1330628
2 chr1  2017424 2017750
3 chr1  2017424 2017750
4 chr1  2018462 2018871
5 chr1  3899411 3899468
6 chr1  4653431 4653724

The expected outcome would be,
head(G3_dat[1:4])

    # A tibble: 6 × 4
      chrom   start     end     dis
      <chr>   <int>   <int>     <int>
    1 chr1  1330426 1330628
    2 chr1  2017424 2017750     686796
    3 chr1  2017424 2017750       -326
    4 chr1  2018462 2018871        712
    5 chr1  3899411 3899468    1880540
    6 chr1  4653431 4653724     753963

There are two functions in the tidyverse suite (specifically dplyr) called lead and lag that allow calculations like this to be done:
library(tidyverse)
G3_dat %>%
   group_by(chrom) %>%
   arrange(start) %>%
   mutate(dis = start - lag(end))

I'm not sure if I've got the lead/lag function correct there (I've updated it based on a suggestion from @ram-rs), but it doesn't take too long to play round with all the combinations to get the right one.
FWIW, %>% (or |> now in native R) is a "pipe" operator that means "take the result of the last function and feed it into the first argument of the next function". It is used extensively in tidyverse scripts. More details here.


text:
Increasing memory usage in AWS instance
I have 16 CPUs on my AWS EC2 instance and I am trying to use kSNP4 to generate an SNP matrix. I was able to successfully run this with over 150 samples - while it took a while to run it did complete - and now I want to run it with over 170 samples but when I use the same command it dies at the second step of the process. I have reason to believe it is due to the processors on my server that I am using more than the tool itself.
Does anyone know the best way to ensure maximum processing ability to run tools on an EC2 instance? Or what I can do to increase the power of the instance?
I know this is more than a technical question than a specific bioinformatic question but I would appreciate any insight!

Do you want to increase CPUs or memory? They are very different things. In either case though, this isn't really something we can help with and you should be talking to your sysadmins or AWS help. Sounds like you need a machine with different specs.

AWS docs outlines how to pick a different instance type for scaling: https://docs.aws.amazon.com/AWSEC2/latest/UserGuide/ec2-instance-resize.html I don't see why a tool that can process 150 samples will crash with 170 unless the complexity goes up geometrically or exponentially somehow

To be honest ... it looks like a bug in the data. 150 -> 170 is not huge and "crashing" is ominous.

You need to diagnose for low memory and consider switching to an R5 processor that prioritises memory. Basically it's easily solvable by switching hardware. There are certainly alot of other approaches but this is specialist HPC.
Personally, I would use Slurm ParallelCluster in AWS. This is very easy to set up from from the configuration file.
The overview of instance types here.
https://aws.amazon.com/ec2/instance-types/


text:
Nextflow: Properly chaining process outputs
I have written a pipeline that is composed of three processes. I seem to be misunderstanding how Nextflow passes the outputs from one process to the next or creates a queue of files. I am using a tumor and paired normal sample from SRA (SRX21185781) to prototype this pipeline.
My first process, trim_QC works fine and outputs the expected number of trimmed fastqs. However, the next process, bwa_align only creates a single bam from one of the two output from the preceding process. My terminal looks like this:
N E X T F L O W  ~  version 23.10.1
Launching `../main.nf` [sick_church] DSL2 - revision: beb03cc0c6

[6e/48d4f6] process > trim_QC (SRR25452858_normal)                  [100%] 2 of 2, cached: 2 ✔
[39/1b8cb2] process > bwa_align (SRR25452858_normal)                [100%] 1 of 1, cached: 1 ✔
[5f/e598c1] process > mapping_stats (SRR25452858_normal_sorted.bam) [100%] 1 of 1, cached: 1 ✔
...

and when running with -with-trace it is even more obvious only one bam file is being processed by the downstream processes:
task_id hash    native_id   name    status  exit    submit  duration    realtime    %cpu    peak_rss    peak_vmem   rchar   wchar
2   d3/0bd651   71127   trim_QC (SRR25452863_tumor) CACHED  0   2024-02-07 16:28:37.449 1h 21m 10s  1h 21m 10s  -   -   -   -   -
1   6e/48d4f6   71128   trim_QC (SRR25452858_normal)    CACHED  0   2024-02-07 16:28:37.459 1h 16m 32s  1h 16m 32s  -   -   -   -   -
3   39/1b8cb2   2754    bwa_align (SRR25452858_normal)  CACHED  0   2024-02-08 10:08:46.033 4h 18m 10s  4h 18m 10s  -   -   -   -   -
4   5f/e598c1   9282    mapping_stats (SRR25452858_normal_sorted.bam)   CACHED  0   2024-02-08 14:26:55.823 53.9s   53.8s   -   -   -   -   -

I am expecting for the bwa-align and mapping_stats processes to also be n of 2.
I have tried explicitly designating the trim_QC.out as a channel (Channel.of(trim_QC.out)) or using an operator to return a queue channel object (i.e., .collect, .collectFile(), etc.) but this either returns errors or fails.
If I try to visualize the outputs from the trim_QC process by adding this to my workflow: trim_QC.out.trimmed_reads.view{it}, it does show both file pairs in my terminal, but only the first one is processed by bwa_align and subsequent steps. Like so:
...
[6e/48d4f6] process > trim_QC (SRR25452858_normal)            [100%] 2 of 2, cached: 2 ✔
[39/1b8cb2] process > bwa_align (SRR25452858_normal)          [100%] 1 of 1, cached: 1 ✔
[5f/e598c1] process > mapping_stats (SRR25452858_normal_so... [100%] 1 of 1, cached: 1 ✔
hello

[SRR25452858_normal, /Users/me/Desktop/fastq2/work/6e/48d4f61cae83ff903abceaadf715e5/SRR25452858_normal_1_val_1.fq.gz, /Users/me/Desktop/fastq2/work/6e/48d4f61cae83ff903abceaadf715e5/SRR25452858_normal_2_val_2.fq.gz]
Analysis complete for SRR25452858_normal_1_val_1.fq.gz
Analysis complete for SRR25452858_normal_2_val_2.fq.gz

[SRR25452863_tumor, /Users/me/Desktop/fastq2/work/d3/0bd6515fcedfc060596b75dd90c28a/SRR25452863_tumor_1_val_1.fq.gz, /Users/me/Desktop/fastq2/work/d3/0bd6515fcedfc060596b75dd90c28a/SRR25452863_tumor_2_val_2.fq.gz]
Analysis complete for SRR25452863_tumor_1_val_1.fq.gz
Analysis complete for SRR25452863_tumor_2_val_2.fq.gz

I am unsure if I am misunderstanding what Nextflow is capable of or how channels work. What I am trying to accomplish is trim adapters and QC result, align trimmed reads followed by indexing, and finally calculate mappings statistics. Any advice is appreciated.
Here is my entire script:
#!/usr/bin/env nextflow
nextflow.enable.dsl = 2

params.fastq = "$projectDir/fastq2/*_{1,2}*"
params.index_dir = "$projectDir/RefGenome/BWAIndex"
params.ref = "hg38.fa"
params.out_dir = "$projectDir/Results"

workflow {

    fastq_ch = Channel.fromFilePairs(params.fastq)
    index_ch = Channel.fromPath(params.index_dir)
    ref_ch = Channel.of(params.ref)

    trim_QC(fastq_ch)
    bwa_align(index_ch, ref_ch, trim_QC.out.trimmed_reads)
    mapping_stats(bwa_align.out.sorted_bam)
}

process trim_QC {
    debug true
    tag "$sample_id"

    publishDir = [
                    [
                        path: { "${params.out_dir}/trimmed/" },
                        mode: 'copy',
                        overwrite: false
                    ],

                    [
                        path: { "${params.out_dir}/stats/" },
                        mode: 'copy',
                        pattern: "*fastqc*",
                        overwrite: false
                    ]
                ]

    input:
    tuple val(sample_id), path(fastq)

    output:
    tuple val("${sample_id}"), path("*_val_1.fq.gz"), path("*_val_2.fq.gz"), emit: trimmed_reads
    path("*fastqc*")

    script:
    def (read1, read2) = fastq
    """
    trim_galore --paired --fastqc ${read1} ${read2}
    """
}

process bwa_align {
    debug true
    tag "$sample_id"

    publishDir "${params.out_dir}/bwa_aligned/", mode: 'copy', overwrite: false

    input:
    path index_dir
    val ref
    tuple val(sample_id), path(read1), path(read2)

    output:
    path("${sample_id}_sorted.bam"), emit: sorted_bam
path("${sample_id}_sorted.bam.bai"), emit: bam_index

    script:
    """
    bwa mem -t 2 ${index_dir}/${ref} ${read1} ${read2} | samtools sort -@2 -o "${sample_id}_sorted.bam" -
samtools index ${sample_id}_sorted.bam > "${sample_id}_sorted.bam.bai"
    """
}

process mapping_stats {
    debug true
    tag "$sorted_bam"

    publishDir "${params.out_dir}/stats/", mode: 'copy', overwrite: false

    input:
    path(sorted_bam)

    output:
    path("${sorted_bam}_mapping_stats.txt")

    script:
    """
    samtools flagstat ${sorted_bam} > "${sorted_bam}_mapping_stats.txt"
    """
}

This is a result of nextflow channel types, queue- and value-channels. Nextflow will make fifo combinations of your channels to run bwa_align until either or both channel are empty. index_ch contains only one item, so once it is used there are no more valid combinations to start new  bwa_align processes.
Either create your index_ch as value channel, or keep the index out of a channel entirely and refer to it as a parameter.

Thanks for your response! I did as you suggested and used Channel.value(params...) for both my index and ref channels. While this did work I still seem to be missing something. While my bwa_align process is now aligning both sets of files, the mapping_stats still is 1 of 1. I was under the impression that any output from a process is a value channel, so this should work. Even more strange is this process only takes a single input. Do you have any advice for how to rectify this or what I am misunderstanding?

Afaik a process will only emit a value channel if ALL its input channels were value channels, otherwise it will be a queue chanel. Make a new question if you have a new question, it is hard to debug your new issue without the exact code. The output you show seems to suggest you are using the -resume option, maybe this is messing with the correct execution after you changed the flow.

Apologies, I was too impatient and responded prematurely. The executor updated to 2 of 2 for the mapping_stats following the completion of the second bwa_align process. Looks like everything is working as expected. Thanks again for your help!


text:
Is ./1 a valid GT value to add to VCF?
I have a VCF file coming from a structural variant detection tool that doesn't add GT value. I can add GT values, for example by using this tool: https://vatools.readthedocs.io/en/latest/vcf_genotype_annotator.html. I know there exist the alternate allele (1), but I don't know if it's homozygous or heterozygous.
My collaborator wants a VCF with genotype. If an SV is called but not genotyped, I was wondering if ./1 could mean unknown allele and alternate allele for a diploid genome (working with human DNA). Although, if it is not a standard format, it may be a problem.
Is ./1 a valid value to add or does it have to be 0/1 or 1/1? Thanks.

No, ./1 isn't valid. You can find the VCF specifications at https://samtools.github.io/hts-specs/. The current version is VCFv4.4 and this is the genotype section which hasn't really changed much from previous versions as far as I recall:

GT (String): Genotype, encoded as allele value preceded by either of /
or | depending on whether that allele is considered phased. The first
phasing indicator may be omitted and is implicitly defined as / if any
phasing indicators are / and | otherwise. The allele values are 0 for
the reference allele (what is in the REF field), 1 for the first
allele listed in ALT, 2 for the second allele list in ALT and so on.
For diploid calls examples could be 0/1, 1 | 0, /0/1, or 1/2, etc.
Haploid calls, e.g. on Y, male non-pseudoautosomal X, or mitochondria,
should be indicated by having only one allele value. A triploid call
might look like 0/0/1, and a partially phased triploid call could be
|0/1/2 to indicate that the first allele is phased with another
variant in the VCF. If a call cannot be made for a sample at a given
locus, ‘.’ must be specified for each missing allele in the GT field
(for example ‘./.’ for a diploid genotype and ‘.’ for haploid
genotype). The meanings of the phasing indicators are as follows (see
the PS and PSL fields below for more details on incorporating phasing
information into the genotypes):

/ : allele is unphased
| : allele is phased (according to the phase-set indicated in PS or PSL)

For symbolic structural variant alleles, GT=0 indicates the absence of any of the ALT symbolic structural
variants defined in the record. Implementer should note that merging a VCF record containing only symbolic
structural variant ALT alleles with a record containing other alleles will result a change of the meaning of the
GT=0 haplotypes from the record containing only symbolic SVs.

This is the bit you are most interested in:

If a call cannot be made for a sample at a given locus, ‘.’ must be specified for each missing allele in the GT field (for example ‘./.’ for a diploid genotype and ‘.’ for haploid genotype).

So, ./1 doesn't mean anything. You could use ./. which is the standard for "not genotyped", or you could simply not have a GT field at all since that isn't actually required by the standard, so a VCF file with no genotype field at all can be a perfectly valid and standard compliant VCF.
In other words, just think about what you really need. Do you absolutely need a GT field? If not, just leave it out and you're done. If you do, you can use ./. if you want to say "I don't know what the genotype here is" (but be aware that some tools interpret ./. as "this might not be a real call"), but don't use ./1 since that won't be understood by anyone.

./1 is valid - it's contained a lot in the output from decompose operations.

@RamRS I don't see why that makes it valid. It isn't mentioned in the specs as far as I can tell, and I don't see anything in the specs that would support that, so why do you say it is valid? There are loads of tools, including commonly used popular ones, that produce invalid vcf files. To be fair, it is an extremely complex format, and it isn't easy to get it right, but the fact that one tool does something is really no indication that that thing is valid. Do you know what they use ./1 to mean?

I assumed that it does not break the format even if it does not adhere to specs 100%, so it's not really invalid. I think they mean that it's pseudo-haploid. When a multiallelic 1/2 GT is broken into 2 biallelic sites, it does not really yield 0/1 and 0/2 but ./1 and ./2 as each entry is kinda sorta haploid.

@RamRS as someone whose day job involves convincing clients that just because their file has a vcf extension, that don't make it no valid vcf, the standard is all I can rely on. If it doesn't comply with the standard, it isn't a valid vcf file and whether or not software can read it correctly is anyone's guess. That said, bcftools, which is famously strict, doesn't complain about 1/. but other tools might. And since it isn't in the specs, we can only say it isn't valid.


text:
GSEA results plot interpretation
I need to plot the results of a GSEA analysis we performed on scRNA data. I am stumbling upon various possibilities for this and sometimes I am struggling a bit to understand how to interpret it. My source is this: https://yulab-smu.top/biomedical-knowledge-mining-book/enrichplot.html
Specifically, in the DotPlot the x axis is geneRatio, in my results table there is no metric called that way, how is it computed? Another question is regarding the ridge plot, the way it is defined is: "The ridgeplot will visualize expression distributions of core enriched genes" This is a bit oscure to me, and in the plot that I generated from my data following that tutoria the x axis is automatically labelled as "enrichment distribution" which in my understanding has a completely different meaning in this case.
Ihave an object of this class: https://rdrr.io/bioc/DOSE/man/gseaResult-class.html, that now I just fed to the functions as explained in that tutorial. The fields of the result table are:
ID     Description setSize enrichmentScore NES pvalue  p.adjust    qvalue  rank    leading_edge    core_enrichment
I don't find any info on how these results are used, I don't know if showing values or plots may be useful as I have more of a methodological doubt.
My question is, if I needed to build these plots on my own without relying on the library documented in the website, starting from the output of cluster profiler, how would I do it? I don't need R code but an understanding of what is done.
Thank you very much in advance.

Hi, it would help if we can see the results table you have. The geneRatio is a measure calculated by clusterProfiler about the number of genes of the gene set found vs those in the gene set.


text:
Oligomer status prediction with Alphafold 2
I am new to using Alphafold 2 and have been primarily using ChimeraX and Colabfold to run some preliminary predictions.
I am particularly interested in using AF2 to predict the homo-oligomerization number of a protein. I have attempted this initially by manually running several predictions with different multiples of the protein end to end ("ABCD:ABCD"). This works but does not provide any direct comparison of plausibility between different oligomerization possibilities, and does not establish a firm prediction as to what the oligomerization state of the protein is.
I have read a paper (https://doi.org/10.1038/s41594-022-00849-w) which appears to use Alphafold2 to make predictions about oligomer count. These researchers appear to use Colabfold for at least part of this calculation.
Can Alphafold 2 (and Colabfold in particular) be used to make predictions about homo-oligomeric state in this way? If so, how?
Thank you.


text:
Different ClinVar classifications/SNPs are linked to the same rsid
I am using R to extract ClinVar classification of SNPs using rsid. Here's an example of my workflow:
Patient data:

Sample
Func
Gene
Codon_Change
dbSNPv151_GRCh38
Chr
Start
End
Ref
Obs

Mary
missense_variant
SPINK5
Cga/Gga
rs199757347
chr5
148091214
148091214
C
G

With this rsid rs199757347, two ClinVar classifications were found with different codon changes, one is C > G (Uncertain significance, missense variant) and the other one C > T (Pathogenic, nonsense variant).
With the rentrez package in R, I can extract both ClinVar classifications:
library(XML)
library(rentrez)

extract_from_esummary(
  entrez_summary(db = "clinvar",
                 id = entrez_search(db = "clinvar", term = "rs199757347")$ids,
                 rettype='xml',
                 parsed = F),
  "germline_classification"
  )

               968176
description    "Uncertain significance"
last_evaluated "2022/01/01 00:00"
review_status  "criteria provided, single submitter"
trait_set      data.frame,2
               351517
description    "Pathogenic"
last_evaluated "2019/01/01 00:00"
review_status  "criteria provided, multiple submitters, no conflicts"
trait_set      data.frame,2

However, only the classification "Uncertain significance" is relevant to my sample, how can I incorporate the alternative alleles in querying ClinVar with rsid? Or do you have a better way to query ClinVar with my sample data?

Do you absolutely need to use rsIDs? The problem isn't clinvar here, it is that the rsID describes two variants.

@terdon, Yes I understand, that's why I wish to "incorporate" the info of ALT allele in the search. And based on my sample data, I think the most useful info for querying ClinVar is the rsID. Do you have other recommendations on querying ClinVar if I don't use rsID? Thank you!

I am afraid not. I am spoiled, the company I work for has a really good API and that's what I use, but it isn't a free product so any answer I post using it would only be helpful to our clients :/ However, my guess is that you should be able to get the nt change if you don't restrict the output to germline_classification only, so then you should be able to parse that to get what you need. Alternatively, if you can convert the rsIDs to HGVS notation (i.e. NM_006846.4(SPINK5):c.652C>G in your case), you can query using that.


text:
After running nf-core, is there a way to map a gene_id to a specific gene's DNA sequence?
I have been running nf-core in Python and it works great! But I have a seemingly simple question that I'm struggling to find an answer for online. After running the nf-core pipeline on my RNA sequencing data, Salmon is able to provide gene counts in a file named salmon.merged.gene_counts.tsv. This file has essentially all the information I am seeking with all counts associated to a unique gene_id. But, for a large number of these gene_ids where no specific gene_name is known, I would like to further find the corresponding DNA sequence (or corresponding transcript's RNA sequence) for each gene_id. Accordingly, my question is: where can I find the information to map the gene_id to the corresponding nucleotides of the gene's (or original transcript's) sequence?
(Apologies in advance if my question is missing information or overly simply and I'm overlooking something obvious — I am quite new to this all.)

Which species? Also:can you provide any examples of gene_ids you have problems with?

Anas platyrhynchos (mallard) — I'm trying to actually do the inverse problem and check if TP53 is present as one of the gene_ids. Looking at related data from https://midb.pnb.uconn.edu/, this gene seems to be missing there as well

re TP53: if one search for p53 gene here: https://www.ensembl.org/Anas_platyrhynchos/Info/Index I do not see the TP53 itself, just the genes binding to it, or downstream targets. If the genome is incomplete/full of gaps then tp53 may be missing from annotation. I guess not just in ENSEMBL but from other genome assembly based as well.

If you look at the orthologues of human TP53 (ENSEMBL), you can find them in chicken, Kakapo, Zebra finch plus a bunch of reptiles. While this is not an ultimate proof that TP53 does exist in your duck, I would expect that some duck RNA-Seq may contain at least fragments of TP53.

Thank you for sharing these thoughts, @darked89. I suppose my simple question is: what is the best way to check the duck RNAseq data for any fragments of TP53? Any specific insights on methods and files to utilize would be helpful!

Additionally, out of curiosity, how is it that the TP53 gene is available on NCBI but unavailable in Ensembl's automated curation pipeline?

re ENSEMBL mallard tp53: from the NCBI there is a link to ENSEMBL rapid. Could be there are two mallard assemblies or at least two versions?

fishing out tp53 RNA-Seq reads

go to NCBI Taxonomy Browser
put Anas platyrhynchos as a search term
click on Anas platyrhynchos
on the right you will have a table with SRA Experiments, click on 11,679
restrict to RNA and paired on the left
send it to Run Selector
order it by the Bases
select few top runs, like SRR16082819

Now it depends what is your exact goal. You may (separate approaches):

get mallard transcriptome as complete as possible, use kallisto just to check that a given pair of FASTQs contains bunch of reads derived from tp53
use STAR to align it to the genome with the contig containing TP53 gene
one can also try to fish out the reads out of FASTQ using tp53 cDNA sequence, but this will likely miss i.e. putative alternative 5- & 3-prime part of the gene

Thank you for this insightful response! For the last part, could you add more details about this step of alignment against the genome? Namely, how exactly can I confirm the contig containing the TP53 gene is present and identified using STAR?

It is this contig NW_024009712 you will have to get both the genome fasta ans GTF annotation from https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_015476345.1/

Also you may try to align the above contig to probably more complete chicken genome to get an idea about how good is the annotation of the p53 gene in mallard/is there anything missing.

Thank you! For clarity, how did you find that contig? Also, strangely the genome fasta file seems to be corrupted in some way (but other files work, e.g. GTF annotation). Is this genome fasta file hosted anywhere else by chance?

I did find this updated version (from 2023-12-14) of the genome here, but is there any way to confirm that the GFF/GTF file includes TP53?

download the annotation then:  zcat GCF_015476345.1_ZJU1.0_genomic.gff.gz  | grep "gene=TP53;" | more

That works great — thank you! Going to test this all out now

I tried converting the GFF to GTF using this method but keep hitting into: KeyError: 'id-LOC110351722' (perhaps due to "gene_id" key not being present for key-value pairings of the .gtf file) — trying to resolve it, but if you have any thoughts, certainly glad to receive it

re conversion gff to gtf: this is a separate question. BTW, why you need it?

It's part of the nf-core pipeline — it expects annotations in GTF format, but the above link only provides GFF files. I tried using gff3_ID_generator.py to generate IDs if the gff file does not have them for every feature, but the error persists

I think I've found the solution here where the larger repo of files are available

Ok. Still, file another question about the gff to gtf conversion of that particular gff file. And it will be helpful if you up vote and accept the answer as a correct one

Yep, just running and testing the pipeline with the updated fasta and gtf files to determine if TP53 is ultimately found with the nf-core pipeline (which I expect it should) — but these new files keep resulting in the pipeline terminating due to an error somewhere in the files which I'm trying to trace down

Got it working now. These gtf files have a lot information that needs to be cleaned / corrected (e.g. missing gene IDs), but it appears to now run!


text:
NEO analysis tutorial and data
I was reading Jason Aten's paper Using genetic markers to orient the edges in quantitative trait networks: The NEO software, which states

The NEO software can be used to orient the edges of gene co-expression networks or quantitative trait networks if the edges can be anchored to genetic marker data. R software tutorials, data, and supplementary material can be downloaded from: http://www.genetics.ucla.edu/labs/horvath/aten/NEO

But it seems the site is down, does anyone know where can I access these data?

Try emailing the author. The horvath UCLA URL is one of the most notorious dead websites.

Like mentioned in the comment above, I would suggest emailing the corresponding author. Scientific articles usually list the email of the corresponding author and this would be the most logical first step in addressing the missing materials.


text:
Why my FeatureScatter looks so different, and what does that mean?
I' am doing a sc-Seq analysis on a dataset from 10x genomics, using Seurat. I followed the standard workflow and I also did it with the SCTransform workflow. The results are very different when I plot a Feature Scatter with the same 2 genes.
Standard workflow:

SCTransform workflow:

So I was wondering why they were so different. I know that the perfect lines in the first one happen because of the normalization and that it usually means that I need more data. But I don't understand why the SCTransform one. looks like it has less data and why it looks so perfectly arranged. My goal is to make a plot that looks normal and natural, and that may or may not have a visual correlation; instead of this weirdly looking and unnatural plots. If you also have a tip for this, I would be very grateful.  Thanks in advance.
Here is my code for the normal workflow:
library(Seurat)
library(tidyverse)
library(ggplot2)
library(gridExtra)
library(harmony)
library(future.apply)
library(cowplot)
library(patchwork)
library("DESeq2")
library(SeuratData)
library(SeuratDisk)
library(SeuratWrappers)
library(sctransform)

data <- Read10X(data.dir = "/Users/rodrigohermoza/Desktop/UTEC/2024-0/PoliSia/rstudio/diferential_expression_2/DE2/filtered_gene_bc_matrices/raw_feature_bc_matrix/")

pbmc <- CreateSeuratObject(counts = data, project = "DE3")
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & percent.mt <15 & nCount_RNA > 50)

pbmc <- NormalizeData(pbmc)
pbmc <- FindVariableFeatures(pbmc, nfeatures = 3000)
pbmc <- ScaleData(pbmc, vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc)

pbmc <- FindNeighbors(pbmc, dims = 1:21)
pbmc <- FindClusters(pbmc, resolution = 1)
pbmc <- RunUMAP(pbmc, dims = 1:21)
DimPlot(pbmc, label = TRUE)

sweep_res <- paramSweep(pbmc,PCs = 1:20, sct = FALSE)
sweep_stat <- summarizeSweep(sweep_res, GT = F)
bcmv <- find.pK(sweep_stat)

ggplot(bcmv, aes(pK, BCmetric, group = 1)) +
  geom_point()+
  geom_line()

#0.005

pK<- bcmv %>%
  filter(BCmetric == max(BCmetric)) %>%
  select(pK)
pK<- as.numeric(as.character(pK[[1]]))

annotations <- pbmc@meta.data$seurat_clusters
homotipic <- modelHomotypic(annotations)
exp <- round(0.08*nrow(pbmc@meta.data))
exp.adj <- round(exp*(1-homotipic))

pbmc <- doubletFinder(pbmc,
                          PCs = 1:20,
                          pK = pK,
                          nExp = exp.adj,
                          reuse.pANN = F,
                          sct = F)
DimPlot(pbmc, reduction = "umap", group.by = "DF.classifications_0.25_0.11_818")

table(pbmc@meta.data$DF.classifications_0.25_0.11_818)
pbmc_filtered <- subset(pbmc, subset = DF.classifications_0.25_0.11_818 == "Singlet")

pbmc.markers <- FindAllMarkers(pbmc_filtered)

FeatureScatter(pbmc_filtered, feature1 = "CD3D", feature2 = "CD3E", slot = "data")

And here is my SCTransform workflow
library(Seurat)
library(tidyverse)
library(ggplot2)
library(gridExtra)
library(harmony)
library(future.apply)
library(cowplot)
library(patchwork)
library("DESeq2")
library(SeuratData)
library(SeuratDisk)
library(SeuratWrappers)
library(sctransform)
library(DoubletFinder)

data <- Read10X(data.dir = "/Users/rodrigohermoza/Desktop/UTEC/2024-0/PoliSia/rstudio/diferential_expression_2/DE2/filtered_gene_bc_matrices/raw_feature_bc_matrix/")
pbmc <- CreateSeuratObject(counts = data, project = "DE3")
pbmc <- PercentageFeatureSet(pbmc, pattern = "^MT-", col.name = "percent.mt")
pbmc <- subset(pbmc, subset = nFeature_RNA > 200 & percent.mt <15 & nCount_RNA > 50)

pbmc <- SCTransform(pbmc, vars.to.regress = "percent.mt")
pbmc <- RunPCA(pbmc)
ElbowPlot(pbmc,ndims = 50)

pbmc <- FindNeighbors(pbmc, dims = 1:20)
pbmc <- FindClusters(pbmc, resolution = 1)
pbmc <- RunUMAP(pbmc, dims = 1:20)
DimPlot(pbmc, label = TRUE)

pbmc.markers <- FindAllMarkers(pbmc)

FeatureScatter(pbmc, feature1 = "CD3D", feature2 = "CD3E", slot = "data", smooth = T)

Check in detail what is plotted. The "standard" is probably some continuous type of data, like scaled counts or something like this, and the second one is probably some sort of more transformed integers, which is not as tightly continuous as the "standard" one. Like imagine you do log2 of 1:10, that gives defined breaks, while something like a Z-score gives a much more continuous distribution without "hard breaks".

Thanks for the clarification :)


text:
Salmon Multiple File SCRIPT giving error
Hi I am using Salmon quantitation for multiple fastq paired files using following code
salmon_index="human_salmon_index"
fastq_dir="fastq/"
for dir in "${fastq_dir}"/SRR*
do
r1_file=$(find "$dir" -name "*_1.fastq")
r2_file=$(find "$dir" -name "*_2.fastq")
samp=$(basename "$dir")
echo "Processing sample ${samp}"
salmon quant -i "$salmon_index" -l A -1 "$r1_file" -2 "$r2_file" -p 10 -o "salmon_out/${samp}_quant"
done

and getting this error after it iterates over all files and start initial processing steps.
Processing sample SRR27756943_2.fastq
Version Info: This is the most recent version of salmon.
### salmon (selective-alignment-based) v1.10.0

Exception : [
The following errors were detected with the read files

======================================================

ERROR: file [] does not appear to exist!

Any solution?

Are you sure the find operations return exactly one file each? The error seems to state that there were no files found for some run. Or it could also be that you have a file that matched your glob pattern that's not a dir, so the find operation fails.
Try this:
salmon_index="human_salmon_index"
fastq_dir=fastq
for dir in ${fastq_dir}/SRR*/
do
r1_file=$(find $dir -name "*_1.fastq")
r2_file=$(find $dir -name "*_2.fastq")
samp=$(basename $dir)
echo "Processing sample ${samp}"
echo salmon quant -i "$salmon_index" -l A -1 "$r1_file" -2 "$r2_file" -p 10 -o "salmon_out/${samp}_quant" | tee salmon_commands.cmd
done

Verify the contents of salmon_commands.cmd then use bash or parallel to run the actual commands.

Hi, I tried but getting this error
find: ‘fastq//SRR*/’: No such file or directory
find: ‘fastq//SRR*/’: No such file or directory
Processing sample SRR*
salmon quant -i ref_salmon/ -l A -1  -2  -p 10 -o salmon_out/SRR*_quant

Try without the quotes. I've made the required changes.


text:
Bedtools genomecov gives me coverage files showing zeros at all nucleotide positions
Also posted on biostars
I have bam files that I've obtained through STAR alignment. Here is the first 10 lines of one of my bam files:
GWNJ-1013:671:GW2304233104th:2:1126:15501:17879 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1126:15130:18427 0       chromosome      1       255     11S139M *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGG  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF::FFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:137        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1206:9272:29857  0       chromosome      1       255     28S108M14S      *       0       0       GGCGGTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGAGATCGGAAGAGCG  FFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1  HI:i:1  AS:i:106        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:30734:10692 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1266:31412:33411 0       chromosome      1       255     24S126M *       0       0       GTTTGGCGTGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTC  FFFFFFFFFFFFFFFF:FFFFFFFFFFFF:FFFF::FFFFFFF:FFFFFFF:FFFFFFFFFFFFFFFF:FFFF,FFFFFFF:FFFFFFFFF:FFFFFFFFFFFFFFF:FFF::FFFFFFFF,FFFFFFFF:FFFFFFFFFFF:F,FFFFF  NH:i:1  HI:i:1  AS:i:124        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1330:21938:10708 0       chromosome      1       255     16S134M *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF,FFFFFFFFFFFFFFFFF:F  NH:i:1  HI:i:1  AS:i:132        nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:17074:4726  0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1367:11939:11584 0       chromosome      1       255     11S82M  *       0       0       ATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGG   FFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFF:F:FFFFFFFFFFFFFFFFFFF  NH:i:1   HI:i:1  AS:i:80 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1372:24740:30843 0       chromosome      1       255     16S48M  *       0       0       TGCGTATTTCCATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAA        FFFFFFFFFFFFFF:FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF        NH:i:1  HI:i:1  AS:i:47 nM:i:0
GWNJ-1013:671:GW2304233104th:2:1401:11107:19319 0       chromosome      1       255     5S144M1S        *       0       0       ATCACAACTGCTCCTCGGAAGTCGACCAACAAGTCAGCTATGACTTGGCATAATTTGTGCCGACAAAATGCGCGCAGAGTATAGGGGTGGATTAACCCCTATTCAACTCTTCGGTAGTGATTTCCGACTTCACGCTACAACAGGAATTGA  FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFFFFFFFFFFFFFFF,FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF:FFF  NH:i:1  HI:i:1  AS:i:142        nM:i:0

This is the script I'm using that contains the genomecov -d -pc -ibam command:
BAM_DIR='./sorted_bam_files'

# Loop over each BAM file in the directory
for file in "$BAM_DIR"/*_high_mapq_reads.sorted.bam
do
    echo "Calculating coverage for $file..."
    # Get the basename of the file
    base=$(basename "$file")
    # Calculate coverage
    bedtools genomecov -d -pc -ibam $file > "${base%.sorted.bam}.coverage"
done

This is the output I keep getting (first 10 lines only shown). I've tried this with at least 30 other bam files but the output is always the same.
chromosome      1       0
chromosome      2       0
chromosome      3       0
chromosome      4       0
chromosome      5       0
chromosome      6       0
chromosome      7       0
chromosome      8       0
chromosome      9       0
chromosome      10      0

Any ideas on why this keeps happening? I don't think there's something wrong with my bam files because I was able to use them for DGE analyses and also checked to see if most reads were getting unmapped but they weren't. Any help is appreciated.

Why is your first chromosome called "chromosome", that's odd.

@ATpoint also second column FLAG == 0 and the fifth MAPQ ==255 meaning: "A value 255 indicates that the mapping quality is not available.". Then the CIGARs looking like the read was indeed mapped 11S82M etc.

@ATpoint The data I have is from bacteria, so I'm only working with 1 chromosome

option -pc is for
-pc     Calculate coverage of pair-end fragments.

but your reads are single-end reads

What tells you they're single ends? Is it because you know this would happen with -pc and single end files or was it something in the bam?

"What tells you they're single ends? Is it b" : the sam flags of each alignemt is zero. Paired would have the '1' bit set . see https://broadinstitute.github.io/picard/explain-flags.html

Using that page, I can see that 11S139M, the first one, is "read paired", "read mapped in proper pair" and "mate unmapped". Could you explain a bit more please? Ideally, adding the detail in the answer so the next user who is confused can also find the info?

the second column is the read flag (section 1.4 of https://samtools.github.io/hts-specs/SAMv1.pdf) . it's always 0  in your example.  0 is single end read mapped on plus strand.

Ah thanks, yes, I was looking at the cigar string, duh. It would be great if you could add some more detail to the answer to make all this clear. Especially since I don't see mention of 0 meaning "single end" either in the specs (but perhaps I missed it) or when using the explain-flags page (which just ignores 0).

"I don't see mention of 0 meaning "single end" either in the specs (but perhaps I missed it)" it's a array of bits https://en.wikipedia.org/wiki/Bit_array . The first bit is for "Read Paired". As it's not set, the read is not part of a "Read Paired", hence it's single end.

@Pierre Some of my bam files have a mix of 0 and 16 for my read flags. So to my understanding I have a mix of single-end and paired-end reverse strand reads. Do I include or omit -pc?


text:
Illumina Final Report seemingly giving impossible genotypes; how do I make sense of this file?
I am working with an Illumina-generated final report file from human genotyping; the MEGA chip I believe. It's in a massive .csv file I've been 'exploring' via unix commands on a compute cluster, something I'm fairly new to.
For a particular SNP (an APOE variant, rs429358, ancestral allele T with the C MAF presumed to be ~15% in this population of ~160 genotyped persons), though, I'm seeing seemingly non-sensical reads. (Table of first two sample IDs' results at bottom.)
I've been trying to make sense of this, and I think the R 1-5 indicates 5 reverse reads and F1-15 indicate 15 forward reads. This also makes sense as to why the R reads would be A-G in the name since that would be the opposing strand of the T-C of the forward reads.
But, as you can see, literally everything on both reads is T-T. This is true for every single person in the sample, everyone is T-T, which even if not for the forward/back mismatch, would still be mathematically virtually impossible.
This is my first time working with raw final report files large enough that I can't simply open them and look at them, and so I'm having trouble determining if there's something wrong with the file, or something wrong with how I'm reading it. Is there any explanation that more experienced people know of for what I'm seeing here? (I will turn to the organization that did the genotyping soon, but I wanted to see if I could figure it out on my own or get as much info as possible to help them troubleshoot first!)

SNP Name
Sample ID
Allele 1 - Forward
Allele 2 - Forward

19:45411941-A-G
001
T
T

19:45411941-A-G-R1
001
T
T

19:45411941-A-G-R2
001
T
T

19:45411941-A-G-R3
001
T
T

19:45411941-A-G-R4
001
T
T

19:45411941-A-G-R5
001
T
T

19:45411941-T-C
001
T
T

19:45411941-T-C-F1
001
T
T

19:45411941-T-C-F10
001
T
T

19:45411941-T-C-F11
001
T
T

19:45411941-T-C-F12
001
T
T

19:45411941-T-C-F13
001
T
T

19:45411941-T-C-F15
001
T
T

19:45411941-T-C-F2
001
T
T

19:45411941-T-C-F3
001
T
T

19:45411941-T-C-F5
001
T
T

19:45411941-T-C-F6
001
T
T

19:45411941-T-C-F7
001
T
T

19:45411941-T-C-F8
001
T
T

19:45411941-T-C-F9
001
T
T

19:45411941-A-G
002
T
T

19:45411941-A-G-R1
002
T
T

19:45411941-A-G-R2
002
T
T

19:45411941-A-G-R3
002
T
T

19:45411941-A-G-R4
002
T
T

19:45411941-A-G-R5
002
T
T

19:45411941-T-C
002
T
T

19:45411941-T-C-F1
002
T
T

19:45411941-T-C-F10
002
T
T

19:45411941-T-C-F11
002
T
T

19:45411941-T-C-F12
002
T
T

19:45411941-T-C-F13
002
T
T

19:45411941-T-C-F15
002
T
T

19:45411941-T-C-F2
002
T
T

19:45411941-T-C-F3
002
T
T

19:45411941-T-C-F5
002
T
T

19:45411941-T-C-F6
002
T
T

19:45411941-T-C-F7
002
T
T

19:45411941-T-C-F8
002
T
T

19:45411941-T-C-F9
002
T
T

Isn't this simply a list of forward and reverse reads that do not support the variant? Can you explain where this table is coming from? How was it generated? What is it supposed to show?

@terdon Thank you, essentially, yes, that did end up being the case. After some reaching out, Illumina explained that with this particular genotyping process they were using at the time (this was done in 2017 I believe) they were attempting to circumvent a common issue with two particular SNPs that are difficult to genotype. They explained that the reason these reads are strange is due to their staff trying new strategies to do the read, but that they strategies did not prove to be successful. So indeed, it was basically an unsupported variant. Thanks!

Oh great, and thanks for reporting back! I will close the question as non-reproducible then since you've figured it out.


text:
How to prove correlation between gene expression and functions using omics data in bioinformatics?
I am new to bioinformatics, coming from a background in mathematical biology. I am currently participating in a project where the Principal Investigator (PI), who is a biologist, has tasked me with proving a correlation between the expression and functions of two or more elements. Specifically, one of them is a kinase enzyme and the other is mTOR. The PI has asked me to use omics data or databases for this purpose.
As someone with limited experience in bioinformatics, I am struggling to understand how a biological pathway or correlation can be definitively proven, beyond basic statistical methods. Additionally, the ultimate goal of this research is to develop treatments using substances that our lab has discovered, which have shown promise in killing specific cancer cells in murine and cell lines (although our data on this is limited).
My question goes as followes- is it possible to demonstrate a connection between these elements using only databases and omics methods? Why might a dynamic modeling approach, such as ordinary differential equations (ODEs) not be enough?
I am feeling quite confused and would greatly appreciate any guidance or direction on how to approach this problem effectively.
Thank you.

In general terms data mining is needed, then supervised learning and finally ODE, because you are performing DE against enzyme kinetics and thats doable, but you'll need a lot of biochemical knowledge.
I would first take the 'omics' data sets and perform a PCA analysis on it and see where your mTOR and kinase genes are located. PCA is known as unsupervised learning and can be accessed by Seurat or better using Scanpy packages. That will help establish the overall behaviour of your key enzymes pre- and post- treatment ... against the control.
Basically you want:

Data sets of cancer mice - pre-treatement
Data set of cancer mice under treatment (drugs targeting kinases)
Healthy controls

See what unsupervised learning looks like between the 3 groups and where your enzymes are. The PCA plots then compared for the 3 groups and you want to see the targets you are interested in changing groups, especially the treatment group shifting between the control and the diseased unsupervised plots. This would provide a basis for exploring supervised learning (outlined below).

Supervised learning What you will do is look at gene expression and train it against your phenotype of interest (cancer and recovery from cancer). The training 'target' is the drug the mice are being treated with for each given data set, obviously they'll be a "untreated" and "control" categories. If your model trains correctly >0.85 to 0.9 accuracy ... post-back and I'll explain how to extract the answer. This approach will extract the biochemical pathway for how the cancer treatment works, but requires more work than described here and that provides the framework for the DE.
Again, a lot of serious knowledge of biochemistry is needed to deploy ODE, unless you've got a brilliant a priori - then you could dive straight in. To do that you'd need to know the exact mechanism of how these enzymes are involved in cancer and then parse that data from your data sets (totally different approach to what is described above). The problem with this approach is if someone knows both the maths and genomics (like a reviewer or examiner) and you've made an error, that will that will be a problem. If you've been through an orthodox approach (above), its a lot stronger justification for the ODE, thus difficult to criticise.
Note The key is what data set types comprise 'omics', I'm assuming this is RNA-based (thats normal), but it might not be.


text:
Extracting P-Values for Co-Expressed Genes Within WGCNA Modules
I'm facing a challenge regarding the extraction of co-expressing genes from each module. After identifying modules, I've been using exportNetworkToCytoscape to visualize genes that co-express together. However, I solely seek the p-values of the co-expressed genes in each module, similar to the format below:
fromNode        toNode      corr-value      p-value
ENSG234        ENSG567       0.90           0.001
ENSG987        ENSG456       0.55           0.03

While the exported data provides useful information such as "fromNode," "toNode,"  "weight", it lacks the p-values (more important) and corr-value  I require for my analysis.
I attempted to use tomSimilarityFromExpr with 'corType = "pearson"' to calculate pairwise correlations between genes, but I'm uncertain about the subsequent steps to obtain the desired p-values.
My goal is to export the co-expression genes of each module, as facilitated by exportNetworkToCytoscape, while also having the corresponding p-values for these co-expressed genes within each module. I'm not interested in linking the modules to external variables; rather, I solely seek the p-values of the co-expressed genes in each module.
Even if there is another way to obtain what I want without going through tomSimilarityFromExpr and exportNetworkToCytoscape, this suggestion would be welcome!
Could you kindly offer guidance on how to proceed in obtaining this information?
Thanks in advance!

The simplest way is to calculate the correlation matrix of genes in the module you want and use that as input to exportNetworkToCytoscape. Something like this:
inMod = moduleColors=="blue"  ## replace 'blue' with the color or module label you want
modExpr = datExpr[, inMod];
modCor = cor(modExpr)
cytoData = exportNetworkToCytoscape(modCor, threshold = 0, ...)
## Replace ... above with appropriate arguments

Then add the p-values to the edge data frame returned within cytoData using function corPvalueStudent, for example as
edgeDF = cytoData$edgeData
nSamples = nrow(datExpr);
edgeDF$pValue = corPvalueStudent(edgeDF$weight, nSamples)

You can then manipulate the correlations and p-values further as necessary, e.g. export them into a csv file.

Thank you,  @Peter. I will try your method, and if it works well, I'll come back to let you know. Thank you in advance.

@Mozart simply mark as "accepted" if this post has been helpful. This looks like good network theory BTW.

Thanks Peter, everything seems to work fine!


text:
How to know if FASTQ/BAM is from reference genome (FASTA)?
I'm new to bioinformatics. I have a problem in which I have a FASTA reference genome and lots of reads in FASTQ files. Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome. I also have access to the BAM/BAM.BAI files associated with each FASTQ, although I can't know how they were aligned.
Thanks

Your post has multiple questions but one of them is not really a question so I'll address everything here:
1. "Some of them could be contaminants, so I'd like to filter them out and get only the ones that align to the reference genome"
Just align the reads. Contaminants should not align well and should be discarded owing to poor mapping quality score.
2. How to figure out if the FASTA file was used for the alignment
This is not really a 100% perfect solution but the BAM header should contain the exact command used for the alignment, which should contain the ref FASTA with path. You can also compare contig names and lengths between the FASTA and the BAM @SQ lines.

"Contaminants should not align well" assumption would not quite hold all the time, for example if the data is from a PDX model involving mouse as well as human reads!

That's true. I was about to suggest something like Xenome or BBSplit but thought these might be simple bacterial contaminants that may not require preprocessing.

In general, this answer wouldn't be more correct than @Ram RS' answer but here is how I would approach:

I would start everything from scratch, at the end of the day you will most probably need to document the exact genome version and its associated annotation file. It looks like you do not have this information at present.

I don't know the experimental set-up but if you are worried about "contaminants", you can have an extra alignment step using BBSplit (link), it is capable of mapping to multiple genomes simultaneously.

I would give a try with nfcore/rnaseq, it abstracts away a lot of complex bioinformatics steps so that you can focus on understanding what these tools do instead of dealing with installing and running these tools. For example you can just use the --bbsplit-fasta-list argument to get rid of your contaminants.