text:
MergeBamAlignment error
I doing the alignment of samples following the GATK pipeline, and doing the MergeBamAlignment,like this:
MergeBamAlignment \
  -ALIGNED $path/file.unsorted.bam \
  -UNMAPPED $path/file.unmapped.bam \
  -O $path/file.merged.bam \
  -R $references/GRCh38.primary_assembly.genome.fa \
  -SO coordinate \
  -TMP_DIR $path/tmp/

I got this error:
Exception in thread "main" java.lang.IllegalArgumentException: Do not use this function to merge dictionaries with different sequences in them. Sequences must be in the same order as well. Found [] and [chr1, chr2, chr3, chr4, chr5, chr6, chr7, chr8, chr9, chr10, chr11, chr12, chr13, chr14, chr15, chr16, chr17, chr18, chr19, chr20, chr21, chr22, chrX, chrY, chrM, GL000008.2, GL000009.2, GL000194.1, GL000195.1, GL000205.2, GL000208.1, GL000213.1, GL000214.1, GL000216.2, GL000218.1, GL000219.1, GL000220.1, GL000221.1, GL000224.1, GL000225.1, GL000226.1, KI270302.1, KI270303.1, KI270304.1, KI270305.1, KI270310.1, KI270311.1, KI270312.1, KI270315.1, KI270316.1, KI270317.1, KI270320.1, KI270322.1, KI270329.1, KI270330.1, KI270333.1, KI270334.1, KI270335.1, KI270336.1, KI270337.1, KI270338.1, KI270340.1, KI270362.1, KI270363.1, KI270364.1, KI270366.1, KI270371.1, KI270372.1, KI270373.1, KI270374.1, KI270375.1, KI270376.1, KI270378.1, KI270379.1, KI270381.1, KI270382.1, KI270383.1, KI270384.1, KI270385.1, KI270386.1, KI270387.1, KI270388.1, KI270389.1, KI270390.1, KI270391.1, KI270392.1, KI270393.1, KI270394.1, KI270395.1, KI270396.1, KI270411.1, KI270412.1, KI270414.1, KI270417.1, KI270418.1, KI270419.1, KI270420.1, KI270422.1, KI270423.1, KI270424.1, KI270425.1, KI270429.1, KI270435.1, KI270438.1, KI270442.1, KI270448.1, KI270465.1, KI270466.1, KI270467.1, KI270468.1, KI270507.1, KI270508.1, KI270509.1, KI270510.1, KI270511.1, KI270512.1, KI270515.1, KI270516.1, KI270517.1, KI270518.1, KI270519.1, KI270521.1, KI270522.1, KI270528.1, KI270529.1, KI270530.1, KI270538.1, KI270539.1, KI270544.1, KI270548.1, KI270579.1, KI270580.1, KI270581.1, KI270582.1, KI270583.1, KI270584.1, KI270587.1, KI270588.1, KI270589.1, KI270590.1, KI270591.1, KI270593.1, KI270706.1, KI270707.1, KI270708.1, KI270709.1, KI270710.1, KI270711.1, KI270712.1, KI270713.1, KI270714.1, KI270715.1, KI270716.1, KI270717.1, KI270718.1, KI270719.1, KI270720.1, KI270721.1, KI270722.1, KI270723.1, KI270724.1, KI270725.1, KI270726.1, KI270727.1, KI270728.1, KI270729.1, KI270730.1, KI270731.1, KI270732.1, KI270733.1, KI270734.1, KI270735.1, KI270736.1, KI270737.1, KI270738.1, KI270739.1, KI270740.1, KI270741.1, KI270742.1, KI270743.1, KI270744.1, KI270745.1, KI270746.1, KI270747.1, KI270748.1, KI270749.1, KI270750.1, KI270751.1, KI270752.1, KI270753.1, KI270754.1, KI270755.1, KI270756.1, KI270757.1].
    at htsjdk.samtools.SAMSequenceDictionary.mergeDictionaries(SAMSequenceDictionary.java:340)
    at picard.sam.SamAlignmentMerger.getDictionaryForMergedBam(SamAlignmentMerger.java:197)
    at picard.sam.AbstractAlignmentMerger.mergeAlignment(AbstractAlignmentMerger.java:384)
    at picard.sam.SamAlignmentMerger.mergeAlignment(SamAlignmentMerger.java:181)
    at picard.sam.MergeBamAlignment.doWork(MergeBamAlignment.java:368)
    at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:309)
    at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:103)
    at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:113)

srun: error: compute-12: task 0: Exited with exit code 1

Before I also confirmed that the files that I'm trying to merge, are with the same reference sequences. So I don't understand what could be the problem.

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.

Welcome @RitaSoares. Firstly is this a human genome? Secondly, what have you done unto this point, ie. how did you create the different alignments, especially what were being used?

Not sure if this is related, but do you get the same error after sorting your aligned bam file? Also, just checking that you aligned the reads to the same reference ($references/GRCh38.primary_assembly.genome.fa)

@jared_mamrot yes but i got a different error after sorting the aligned bam files, with SortSamSpark: "org.broadinstitute.hellbender.tools.spark.pipelines. SortSamSpark done. Elapsed time: 0.07 minutes.Runtime.totalMemory()=2854748160
htsjdk.samtools.util.RuntimeEOFException: Premature EOF. Expected 4 but only received 0; BinaryCodec in readmode; streamed file (filename not available)". But according to this error that happend to me before, normally i put in the code of MergeBamAlignment the argument "SO coordinate"and  this corrects my problem but with this samples doesn't happen that.

Sounds to me like your sample might be corrupted/truncated. What does picard ValidateSamFile output for that sample?

I also think that is the case... i will try to do ValidateSamFile to make sure. Thank for your answer


text:
Convert ensemble genes to gene names in my sigle cell signature matrix annotated with seurat in R?
I have single cell data which I want to annotate with the seurat package, so I ran this code in R:
lung.data <- Read10X(data.dir ="C:/Users/Desktop/SINGLE-CELL/DATA",gene.column = 1)

lung <- CreateSeuratObject(counts = lung.data, project = "lung", min.cells = 3, min.features = 200)

And the output of the features of this seurat object is:
    >rownames(lung)
    [1] "ENSG00000000003" "ENSG00000000419"
 [3] "ENSG00000000457" "ENSG00000000460"
 [5] "ENSG00000000938" "ENSG00000000971"
 [7] "ENSG00000001036" "ENSG00000001084"
 [9] "ENSG00000001167" "ENSG00000001460"
[11] "ENSG00000001461" "ENSG00000001497"
[13] "ENSG00000001561" "ENSG00000001617"
[15] "ENSG00000001626" "ENSG00000001629"
[17] "ENSG00000001630" "ENSG00000001631"
[19] "ENSG00000002016" "ENSG00000002079"

How can I change the row names of this Seurat Object from Ensembl to gene name, knowing that not all the ensemble genes have conversion to gene name, in order to obtain the seurat object with the gene names as row names and cell names as column names?

What is a "gene name"? I checked one of your ENSGs and it looks like a human gene, so does that mean you are only working with human sequences? And, if so, do you then mean the HGNC symbol for the gene? Or do you also need to handle other species?

From Seurat developers:

the functions Read10X has an argument gene.column that can be used to change between using the gene name and the ensembl ID when reading in 10x Genomics datasets. You can also use the use.names parameter when using Read10X_h5."

I think BioMart will do the gene ID -> gene name conversion you're after:
https://www.ensembl.org/info/data/biomart/index.html
To use BioMart to convert the gene IDs to gene names, you can use the BioMaRt bioconductor R package. Installation and running instructions can be found here: ensembl.org/info/data/biomart/biomart_r_package.html
For this query, you will need to use the human genes dataset:
ensembl = useEnsembl(biomart="ensembl", dataset="hsapiens_gene_ensembl")

Then use your list of gene IDs as the 'Filter', with list of IDs as the values:
filters = 'ensembl_gene_id'

Most human protein-coding genes have an associated HGNC symbol, and these are used in Ensembl. Therefore, you can select the 'gene_name' attribute to retrieve the human gene names, and this will usually correspond to the HGNC gene names, as stored in the Ensembl database.
More information on gene names can be found in the Ensembl documentation: https://ensembl.org/info/genome/genebuild/gene_names.html

source: Stack Exchange

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

sort: votes


text:
ERROR IN FUNCOTATOR: [SQLITE_IOERR_LOCK] I/O error in the advisory file locking logic (disk I/O error)
I ran the "funcotator" step from the GATK somatic short variant discovery (snvs + indels) pipeline and used the output of filtermutectcalls as input, but I obtain this error:
17:00:18.816 INFO  NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/gatk/gatk-package-4.3.0.0-local.jar!/com/intel/gkl/native/libgkl_compression.so
17:00:19.105 INFO  Funcotator - ------------------------------------------------------------
17:00:19.106 INFO  Funcotator - The Genome Analysis Toolkit (GATK) v4.3.0.0
17:00:19.106 INFO  Funcotator - For support and documentation go to https://software.broadinstitute.org/gatk/
17:00:19.106 INFO  Funcotator - Executing as rita.soares@compute-19 on Linux v5.4.0-148-generic amd64
17:00:19.106 INFO  Funcotator - Java runtime: OpenJDK 64-Bit Server VM v1.8.0_242-8u242-b08-0ubuntu3~18.04-b08
17:00:19.106 INFO  Funcotator - Start Date/Time: November 24, 2023 5:00:18 PM GMT
17:00:19.106 INFO  Funcotator - ------------------------------------------------------------
17:00:19.106 INFO  Funcotator - ------------------------------------------------------------
17:00:19.107 INFO  Funcotator - HTSJDK Version: 3.0.1
17:00:19.107 INFO  Funcotator - Picard Version: 2.27.5
17:00:19.107 INFO  Funcotator - Built for Spark Version: 2.4.5
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.COMPRESSION_LEVEL : 2
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_READ_FOR_SAMTOOLS : false
17:00:19.107 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_SAMTOOLS : true
17:00:19.108 INFO  Funcotator - HTSJDK Defaults.USE_ASYNC_IO_WRITE_FOR_TRIBBLE : false
17:00:19.108 INFO  Funcotator - Deflater: IntelDeflater
17:00:19.108 INFO  Funcotator - Inflater: IntelInflater
17:00:19.108 INFO  Funcotator - GCS max retries/reopens: 20
17:00:19.108 INFO  Funcotator - Requester pays: disabled
17:00:19.108 INFO  Funcotator - Initializing engine
17:00:19.577 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FilterMutectCalls/10LC_tumor_S2.filtered.vcf.gz
17:00:19.880 INFO  Funcotator - Done initializing engine
17:00:19.880 INFO  Funcotator - Skipping sequence dictionary validation.
17:00:19.880 INFO  Funcotator - Processing user transcripts/defaults/overrides...
17:00:19.881 INFO  Funcotator - Initializing data sources...
17:00:19.883 INFO  DataSourceUtils - Initializing data sources from directory: /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s
17:00:19.953 INFO  DataSourceUtils - Data sources version: 1.7.2020429s
17:00:19.953 INFO  DataSourceUtils - Data sources source: ftp://gsapubftp-anonymous@ftp.broadinstitute.org/bundle/funcotator/funcotator_dataSources.v1.7.20200429s.tar.gz
17:00:19.953 INFO  DataSourceUtils - Data sources alternate source: gs://broad-public-datasets/funcotator/funcotator_dataSources.v1.7.20200429s.tar.gz
17:00:20.435 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.annotation.REORDERED.gtf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:20.437 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.pc_transcripts.fa -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.pc_transcripts.fa
17:00:20.513 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:20.584 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Familial_Cancer_Genes.no_dupes.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/familial/hg38/Familial_Cancer_Genes.no_dupes.tsv
17:00:20.711 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/dnaRepairGenes.20180524T145835.csv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dna_repair_genes/hg38/dnaRepairGenes.20180524T145835.csv
17:00:20.779 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_fusion.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_fusion/hg38/cosmic_fusion.tsv
17:00:20.911 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hgnc_download_Nov302017.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/hgnc/hg38/hgnc_download_Nov302017.tsv
17:00:20.967 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/simple_uniprot_Dec012014.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/simple_uniprot/hg38/simple_uniprot_Dec012014.tsv
17:00:21.034 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/CancerGeneCensus_Table_1_full_2012-03-15.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cancer_gene_census/hg38/CancerGeneCensus_Table_1_full_2012-03-15.txt
17:00:21.073 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/achilles_lineage_results.import.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/achilles/hg38/achilles_lineage_results.import.txt
17:00:21.161 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:21.171 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xhgnc_v90_38.hg38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xhgnc/hg38/gencode_xhgnc_v90_38.hg38.tsv
17:00:21.187 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Cosmic.db -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic/hg38/Cosmic.db
17:00:21.341 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:21.357 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xrefseq_v90_38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xrefseq/hg38/gencode_xrefseq_v90_38.tsv
17:00:21.378 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_tissue.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_tissue/hg38/cosmic_tissue.tsv
17:00:21.379 INFO  Funcotator - Finalizing data sources (this step can be long if data sources are cloud-based)...
17:00:21.379 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.annotation.REORDERED.gtf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:21.380 INFO  DataSourceUtils - Setting lookahead cache for data source: Gencode : 100000
17:00:21.553 INFO  FeatureManager - Using codec GencodeGtfCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.annotation.REORDERED.gtf
17:00:21.691 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode.v34.pc_transcripts.fa -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode/hg38/gencode.v34.pc_transcripts.fa
17:00:27.190 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.190 INFO  DataSourceUtils - Setting lookahead cache for data source: dbSNP : 100000
17:00:27.412 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.563 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hg38_All_20180418.vcf.gz -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.666 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dbsnp/hg38/hg38_All_20180418.vcf.gz
17:00:27.814 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Familial_Cancer_Genes.no_dupes.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/familial/hg38/Familial_Cancer_Genes.no_dupes.tsv
17:00:27.837 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/dnaRepairGenes.20180524T145835.csv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/dna_repair_genes/hg38/dnaRepairGenes.20180524T145835.csv
17:00:27.853 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_fusion.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_fusion/hg38/cosmic_fusion.tsv
17:00:27.875 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/hgnc_download_Nov302017.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/hgnc/hg38/hgnc_download_Nov302017.tsv
17:00:28.142 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/simple_uniprot_Dec012014.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/simple_uniprot/hg38/simple_uniprot_Dec012014.tsv
17:00:28.400 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/CancerGeneCensus_Table_1_full_2012-03-15.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cancer_gene_census/hg38/CancerGeneCensus_Table_1_full_2012-03-15.txt
17:00:28.426 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/achilles_lineage_results.import.txt -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/achilles/hg38/achilles_lineage_results.import.txt
17:00:28.449 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.449 INFO  DataSourceUtils - Setting lookahead cache for data source: ClinVar_VCF : 100000
17:00:28.502 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.702 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/clinvar_20180429_hg38.vcf -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.780 INFO  FeatureManager - Using codec VCFCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/clinvar/hg38/clinvar_20180429_hg38.vcf
17:00:28.835 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xhgnc_v90_38.hg38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xhgnc/hg38/gencode_xhgnc_v90_38.hg38.tsv
17:00:29.651 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/Cosmic.db -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic/hg38/Cosmic.db
17:00:29.807 INFO  DataSourceUtils - Setting lookahead cache for data source: Oreganno : 100000
17:00:29.828 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:29.871 INFO  FeatureManager - Using codec XsvLocatableTableCodec to read file file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.config
17:00:30.155 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
17:00:30.158 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/oreganno.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/oreganno/hg38/oreganno.tsv
WARNING 2023-11-24 17:00:30     AsciiLineReader Creating an indexable source for an AsciiFeatureCodec using a stream that is neither a PositionalBufferedStream nor a BlockCompressedInputStream
17:00:30.168 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/gencode_xrefseq_v90_38.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/gencode_xrefseq/hg38/gencode_xrefseq_v90_38.tsv
17:00:30.284 INFO  DataSourceUtils - Resolved data source file path: file:///home/rita.soares/cosmic_tissue.tsv -> file:///mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s/cosmic_tissue/hg38/cosmic_tissue.tsv
17:00:30.398 INFO  Funcotator - Initializing Funcotator Engine...
17:00:30.400 INFO  FuncotatorEngine - Using given VCF and Reference.  No conversion required.
17:00:30.402 INFO  Funcotator - Creating a MAF file for output: file:/mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FUNCOTATOR/10LC_tumor_S2.FiltPASS.maf
17:00:30.420 INFO  ProgressMeter - Starting traversal
17:00:30.420 INFO  ProgressMeter -        Current Locus  Elapsed Minutes    Variants Processed  Variants/Minute
17:00:42.145 INFO  ProgressMeter -         chr1:1129652              0.2                  1000           5117.3
17:00:52.621 INFO  ProgressMeter -         chr1:1544800              0.4                  2000           5405.4
17:01:03.582 INFO  ProgressMeter -         chr1:2088265              0.6                  3000           5427.9
17:01:23.573 INFO  ProgressMeter -         chr1:2868977              0.9                  5000           5644.2
17:01:45.697 INFO  ProgressMeter -         chr1:3708722              1.3                  7000           5579.4
17:01:47.516 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000581512.2 for variant: chr1:3783110-3783112(TAA* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000581512.2  Variant: [VC Unknown @ chr1:3783110-3783112 Q. of type=INDEL alleles=[TAA*, T] attr={AS_FilterStatus=[weak_evidence, contamination], AS_SB_TABLE=[0, 0|0, 0], DP=1, ECNT=2, GERMQ=1, MBQ=[0, 31], MFRL=[0, 259], MMQ=[60, 60], MPOS=66, POPAF=2.21, ROQ=93, RPA=[15, 13], RU=A, STR=true, STRQ=93, TLOD=4.20} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:PGT:PID:PS:SB     0|1:0,1:0.667:1:0,0:0,1:0,1:0|1:3783099_C_T:3783099:0,0,1,0 filters=haplotype,germline,contamination,weak_evidence

... (CODE OMITTED BECAUSE OF CHARACTERS LIMINATION OF THE BODY OF THE QUESTION)

17:06:03.884 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000375071.4 for variant: chr1:20589096-20589097(TC* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000375071.4  Variant: [VC Unknown @ chr1:20589096-20589097 Q. of type=INDEL alleles=[TC*, T] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[108, 95|79, 72], DP=369, ECNT=2, GERMQ=1, MBQ=[20, 20], MFRL=[119, 123], MMQ=[60, 60], MPOS=26, PON=true, POPAF=0.256, ROQ=93, RPA=[3, 2], RU=C, STR=true, STRQ=93, TLOD=206.22} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB   0/1:203,151:0.422:354:57,38:55,42:114,83:108,95,79,72 filters=panel_of_normals,germline
17:06:03.885 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000375071.4 for problem variant: chr1:20589096-20589097(TC* -> T)
17:06:20.587 INFO  ProgressMeter -        chr1:21789645              5.8                 38000           6511.2
17:06:38.318 INFO  ProgressMeter -        chr1:22910935              6.1                 40000           6523.5
17:06:55.667 INFO  ProgressMeter -        chr1:24196391              6.4                 42000           6541.3
17:07:12.288 INFO  ProgressMeter -        chr1:25333744              6.7                 44000           6569.3
17:07:29.628 INFO  ProgressMeter -        chr1:26561968              7.0                 46000           6583.9
17:08:54.624 INFO  ProgressMeter -        chr1:27258140              8.4                 47000           5593.0
17:09:11.523 INFO  ProgressMeter -        chr1:28509457              8.7                 49000           5641.9
17:09:14.750 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000294409.2 for variant: chr1:28714873-28714874(TA* -> T): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000294409.2  Variant: [VC Unknown @ chr1:28714873-28714874 Q. of type=INDEL alleles=[TA*, T] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[1, 8|0, 5], DP=14, ECNT=1, GERMQ=1, MBQ=[29, 29], MFRL=[277, 190], MMQ=[60, 60], MPOS=19, PON=true, POPAF=1.96, ROQ=93, RPA=[13, 12], RU=A, STR=true, STRQ=1, TLOD=6.90} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB   0/1:9,5:0.377:14:1,1:4,1:8,5:1,8,0,5 filters=panel_of_normals,germline,slippage
17:09:14.750 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000294409.2 for problem variant: chr1:28714873-28714874(TA* -> T)
17:09:28.156 INFO  ProgressMeter -        chr1:29742106              9.0                 51000           5690.5
17:09:44.363 INFO  ProgressMeter -        chr1:31000060              9.2                 53000           5740.7
17:10:02.660 INFO  ProgressMeter -        chr1:32309371              9.5                 55000           5766.8
17:10:21.404 INFO  ProgressMeter -        chr1:33875394              9.8                 57000           5787.0
17:10:39.604 INFO  ProgressMeter -        chr1:35564548             10.2                 59000           5811.1
17:11:13.558 INFO  ProgressMeter -        chr1:36959618             10.7                 61000           5690.8
17:11:30.930 INFO  ProgressMeter -        chr1:38301837             11.0                 63000           5722.9
17:11:48.409 INFO  ProgressMeter -        chr1:39672835             11.3                 65000           5752.3
17:12:04.997 INFO  ProgressMeter -        chr1:40888483             11.6                 67000           5787.7
17:12:22.171 INFO  ProgressMeter -        chr1:42389559             11.9                 69000           5816.6
17:13:07.894 INFO  ProgressMeter -        chr1:43856329             12.6                 71000           5624.0
17:13:25.958 INFO  ProgressMeter -        chr1:45154774             12.9                 73000           5647.7
17:13:43.676 INFO  ProgressMeter -        chr1:46612646             13.2                 75000           5672.8
17:14:00.829 INFO  ProgressMeter -        chr1:48197446             13.5                 77000           5700.8
17:14:19.050 INFO  ProgressMeter -        chr1:50725644             13.8                 79000           5720.3
17:14:29.607 ERROR GencodeFuncotationFactory - Problem creating a GencodeFuncotation on transcript ENST00000361556.9 for variant: chr1:51729808-51729817(GGGGGGTGAA* -> G): Variant overlaps transcript but is not completely contained within it. Funcotator cannot currently handle this case. Transcript: ENST00000361556.9  Variant: [VC Unknown @ chr1:51729808-51729817 Q. of type=INDEL alleles=[GGGGGGTGAA*, G] attr={AS_FilterStatus=SITE, AS_SB_TABLE=[44, 1|46, 0], DP=110, ECNT=1, GERMQ=1, MBQ=[37, 37], MFRL=[163, 167], MMQ=[60, 60], MPOS=35, POPAF=1.89, ROQ=93, RPA=[2, 1], RU=GGGGGTGAA, STR=true, STRQ=93, TLOD=175.52} GT=GT:AD:AF:DP:F1R2:F2R1:FAD:SB 0/1:45,46:0.514:91:23,24:12,22:44,46:44,1,46,0 filters=germline
17:14:29.607 WARN  GencodeFuncotationFactory - Creating default GencodeFuncotation on transcript ENST00000361556.9 for problem variant: chr1:51729808-51729817(GGGGGGTGAA* -> G)
17:14:29.715 INFO  ProgressMeter -        chr1:51744810             14.0                 80000           5719.1
17:14:47.317 INFO  ProgressMeter -        chr1:53112621             14.3                 82000           5741.6
17:15:03.760 INFO  ProgressMeter -        chr1:54313159             14.6                 84000           5770.9
17:15:20.710 INFO  ProgressMeter -        chr1:55702719             14.8                 86000           5795.9
17:15:38.363 INFO  ProgressMeter -        chr1:57618601             15.1                 88000           5815.3
17:15:55.381 INFO  ProgressMeter -        chr1:59379903             15.4                 90000           5838.1
17:16:13.770 INFO  ProgressMeter -        chr1:61650616             15.7                 92000           5851.5
17:16:29.978 INFO  ProgressMeter -        chr1:63193764             16.0                 94000           5877.7
17:17:46.249 INFO  ProgressMeter -        chr1:64973681             17.3                 96000           5560.8
17:18:29.297 INFO  ProgressMeter -        chr1:65743799             18.0                 97000           5394.5
17:18:44.743 INFO  VcfFuncotationFactory - ClinVar_VCF 20180429_hg38 cache hits/total: 0/138
17:18:44.746 INFO  VcfFuncotationFactory - dbSNP 9606_b151 cache hits/total: 0/63294
17:18:44.748 INFO  Funcotator - Shutting down engine
[November 24, 2023 5:18:44 PM GMT] org.broadinstitute.hellbender.tools.funcotator.Funcotator done. Elapsed time: 18.43 minutes.
Runtime.totalMemory()=2054684672

org.broadinstitute.hellbender.exceptions.GATKException: Unable to query the database for geneName: LINC01702
        at org.broadinstitute.hellbender.tools.funcotator.dataSources.cosmic.CosmicFuncotationFactory.createFuncotationsOnVariant(CosmicFuncotationFactory.java:320)
        at org.broadinstitute.hellbender.tools.funcotator.DataSourceFuncotationFactory.determineFuncotations(DataSourceFuncotationFactory.java:245)
        at org.broadinstitute.hellbender.tools.funcotator.DataSourceFuncotationFactory.createFuncotations(DataSourceFuncotationFactory.java:211)
        at org.broadinstitute.hellbender.tools.funcotator.FuncotatorEngine.createFuncotationMapForVariant(FuncotatorEngine.java:178)
        at org.broadinstitute.hellbender.tools.funcotator.Funcotator.enqueueAndHandleVariant(Funcotator.java:924)
        at org.broadinstitute.hellbender.tools.funcotator.Funcotator.apply(Funcotator.java:878)
        at org.broadinstitute.hellbender.engine.VariantWalker.lambda$traverse$0(VariantWalker.java:104)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.accept(ForEachOps.java:183)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.stream.ReferencePipeline$2$1.accept(ReferencePipeline.java:175)
        at java.util.stream.ReferencePipeline$3$1.accept(ReferencePipeline.java:193)
        at java.util.Iterator.forEachRemaining(Iterator.java:116)
        at java.util.Spliterators$IteratorSpliterator.forEachRemaining(Spliterators.java:1801)
        at java.util.stream.AbstractPipeline.copyInto(AbstractPipeline.java:482)
        at java.util.stream.AbstractPipeline.wrapAndCopyInto(AbstractPipeline.java:472)
        at java.util.stream.ForEachOps$ForEachOp.evaluateSequential(ForEachOps.java:150)
        at java.util.stream.ForEachOps$ForEachOp$OfRef.evaluateSequential(ForEachOps.java:173)
        at java.util.stream.AbstractPipeline.evaluate(AbstractPipeline.java:234)
        at java.util.stream.ReferencePipeline.forEach(ReferencePipeline.java:485)
        at org.broadinstitute.hellbender.engine.VariantWalker.traverse(VariantWalker.java:102)
        at org.broadinstitute.hellbender.engine.GATKTool.doWork(GATKTool.java:1095)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:140)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
        at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
        at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
        at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
        at org.broadinstitute.hellbender.Main.main(Main.java:289)
Caused by: org.sqlite.SQLiteException: [SQLITE_IOERR_LOCK]  I/O error in the advisory file locking logic (disk I/O error)
        at org.sqlite.core.DB.newSQLException(DB.java:1030)
        at org.sqlite.core.DB.newSQLException(DB.java:1042)
        at org.sqlite.core.DB.execute(DB.java:884)
        at org.sqlite.core.CoreStatement.exec(CoreStatement.java:80)
        at org.sqlite.jdbc3.JDBC3Statement.executeQuery(JDBC3Statement.java:68)
        at org.broadinstitute.hellbender.tools.funcotator.dataSources.cosmic.CosmicFuncotationFactory.createFuncotationsOnVariant(CosmicFuncotationFactory.java:288)
        ... 26 more
Using GATK jar /gatk/gatk-package-4.3.0.0-local.jar
Running:
    java -Dsamjdk.use_async_io_read_samtools=false -Dsamjdk.use_async_io_write_samtools=true -Dsamjdk.use_async_io_write_tribble=false -Dsamjdk.compression_level=2 -jar /gatk/gatk-package-4.3.0.0-local.jar Funcotator -R /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/bwa_prefiles/GRCh38.primary_assembly.genome.fa -V /mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FilterMutectCalls/10LC_tumor_S2.filtered.vcf.gz -O /mnt/nfs/lobo/IMM-NFS/cfranco/MCPP-Lung/DNA/Variant_Calling/FUNCOTATOR/10LC_tumor_S2.FiltPASS.maf --output-file-format MAF --ref-version hg38 --data-sources-path /mnt/nfs/lobo/SALMEIDA-NFS/lcosta/MCPP-Lung/gatk_files/funcotator_dataSources.v1.7.20200521s --disable-sequence-dictionary-validation
srun: error: compute-19: task 0: Exited with exit code 3
srun: launch/slurm: _step_signal: Terminating StepId=2949215.0

Are you sure there's a problem? The error appears after the Funcotator done. Elapsed time: 18.43 minutes. message, which suggests it had actually finished. Then, the error itself (Unable to query the database for gene... and Caused by: [...]  I/O error in the advisory file locking logic (disk I/O error)) suggests a problem with reading some DB that may well be transient. Did you try running it again?

Yes i did runned again but i got the same error and i only obtain a vcf file and not a maf file as output of Funcotator

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


text:
How can I download multiple SRR* files at a time when using sratoolkit's fastq-dump?
I have multiple SRR files to download, and their accession is in a list:
SRR9004325
SRR9004326
SRR9004327

I read the NCBI document but don't find how to achieve this using the command "fastq-dump".
And I noticed that if I use the command prefetch, I can download from the list file:
prefetch --option-file SraAccList.txt

So is there any parameter in "fastq-dump" to download accessions stored in a file?
As @terdon said, it is better to use fasterq-dump, and the command:
fasterq-dump --option-file SraAccList.txt doesn't work. When I try it, I just get the usage message:
$ fasterq-dump --option-file SraAccList.txt

Usage:
  fasterq-dump <path> [options]
  fasterq-dump <accession> [options]

Options:
  -F|--format                      format (special, fastq, default=fastq)
  -o|--outfile                     output-file
  -O|--outdir                      output-dir
  -b|--bufsize                     size of file-buffer dflt=1MB
  -c|--curcache                    size of cursor-cache dflt=10MB
  -m|--mem                         memory limit for sorting dflt=100MB
  -t|--temp                        where to put temp. files dflt=curr dir
  -e|--threads                     how many thread dflt=6
  -p|--progress                    show progress
  -x|--details                     print details
  -s|--split-spot                  split spots into reads
  -S|--split-files                 write reads into different files
  -3|--split-3                     writes single reads in special file
  --concatenate-reads              writes whole spots into one file
  -Z|--stdout                      print output to stdout
  -f|--force                       force to overwrite existing file(s)
  --skip-technical                 skip technical reads
  --include-technical              include technical reads
  -M|--min-read-len                filter by sequence-len
  --table                          which seq-table to use in case of pacbio
  -B|--bases                       filter by bases
  -A|--append                      append to output-file
  --fasta                          produce FASTA output
  --fasta-unsorted                 produce FASTA output, unsorted
  --fasta-ref-tbl                  produce FASTA output from REFERENCE tbl
  --fasta-concat-all               concatenate all rows and produce FASTA
  --internal-ref                   extract only internal REFERENCEs
  --external-ref                   extract only external REFERENCEs
  --ref-name                       extract only these REFERENCEs
  --ref-report                     enumerate references
  --use-name                       print name instead of seq-id
  --seq-defline                    custom defline for sequence:  $ac=accession,
                               $sn=spot-name,  $sg=spot-group, $si=spot-id,
                                   $ri=read-id, $rl=read-length
  --qual-defline                   custom defline for qualities:  same as
                                   seq-defline
  -U|--only-unaligned              process only unaligned reads
  -a|--only-aligned                process only aligned reads
  --disk-limit                     explicitly set disk-limit
  --disk-limit-tmp                 explicitly set disk-limit for temp. files
  --size-check                     switch to control: on=perform size-check
                                   (default),  off=do not perform size-check,
                                   only=perform size-check only
  --ngc <PATH>                     PATH to ngc file

  -h|--help                        Output brief explanation for the program.
  -V|--version                     Display the version of the program then
                                   quit.
  -L|--log-level <level>           Logging level as number or enum string. One
                                   of (fatal|sys|int|err|warn|info|debug) or
                                   (0-6) Current/default is warn.
  -v|--verbose                     Increase the verbosity of the program
                                   status messages. Use multiple times for more
                                   verbosity. Negates quiet.
  -q|--quiet                       Turn off all status messages for the
                                   program. Negated by verbose.
  --option-file <file>             Read more options and parameters from the
                                   file.
for more information visit:
   https://github.com/ncbi/sra-tools/wiki/HowTo:-fasterq-dump
   https://github.com/ncbi/sra-tools/wiki/08.-prefetch-and-fasterq-dump

fasterq-dump : 3.0.6

I don't understand. What is wrong with prefetch? That's exactly what you need, right? fastq-dump (which I believe is obsolete, you should use fasterq-dump instead) is to extract fastq files from the SRR, not for downloading.

@terdon From the documentation: " You can also avoid the prefetch step and download and convert the Run in one step by entering just the Run accession without the .sra extension in your fasterq-dump or sam-dump command".  I just want to download the fastq file.

OK, so what happens if you run fasterq-dump --option-file SraAccList.txt? If you look at the tool's options (run fasterq-dump with no arguments), it shows that the option is there.

Hmm. I tried that and it doesn't seem to work. But please add that to the question. Have you gotten the command to work correctly for a single file? If so, you could just loop over the accessions.

@terdon Yes, the options " --option-file" is there, but the description about it is "Read more options and parameters from the file". Emmm, it is a little confusing.

If fasterq-dump can't read from a file, you can simply loop over the file yourself:
while read srr; do
  fasterq-dump "$srr" -o "$srr".fastq
done < SraAccList.txt

That will iterate over your SraAccList.txt file and extract the fastqs for each accession. I tested with this command:
fasterq-dump SRR9004325 -o ff.fastq

And that created files ff_1.fastq and ff_2.fastq, so you can expect files named like SRR9004325_1.fastq and SRR9004325_1.fastq if you use the loop above.
You could try and parallelize it but I doubt it is worth it since the bottleneck will always be the download speed.

Note that I have switched from getting fastq files to first downloading SRA files and then converting them to fastq myself, as this is (in my experience), much faster (SRA files are smaller and somehow downloading fastq files directly was very slow for me).
I first do this:
prefetch {sample} -O prefetch_dir --max-size 50G
Then
fasterq-dump {sra_file} --outdir fastq_raw_dir
gzip fastq_raw_dir/{sample}_1.fastq
gzip fastq_raw_dir/{sample}_2.fastq

I use Snakemake to parallelize.


text:
How to analyze GC content of fasta file?
I have the assembly file of a certain lizard species, I wanted to recreate Fig. 2 of the following paper, which compares the GC content distribution over 5kb windows, how would I go about doing this?
https://bmcecolevol.biomedcentral.com/articles/10.1186/s12862-022-02086-7#Sec10

I would suggest using bedtools makewindows and bedtools nuc to compute the GC content for genomic windows as shown in Fig 2 of the linked paper.
See related question here.
This will slice the genome into 5Kbp (or whatever size specified to makewindows) pieces and compute the GC content of each piece, which can be plotted as a distribution as in the figure.
Something like this would work for the 5000 bp example:
# assemble a bedtools compatible genome file:
# see https://www.biostars.org/p/70795/ for hints how to do this.
bedtools makewindows -w 5000 -g my_genome.genome > windows.bed # you just made this .genome file manually
bedtools nuc -bed windows.bed -fi my_genome.fasta > windows_gc.bed

The bed output file will contain the GC content of each window in the 7_pct_gc column (according to my understanding of the output file).

Please note that you do not have to make ".genome" file manually if you have samtools. You could simply generate an index file by running samtools faidx genome.fasta, and use the output file with extension .fasta.fai as input in your makewindwos command.

tail -n +2 windows_gc.bed | cut -f 5 | awk '{printf("%d\n",$1*100)}' | sort | uniq -c | sed 's/^ *//'| sed 's/ /,/'>output.csv. This should print a CSV file which you can easily plot.

Try this small package/library:
https://github.com/willros/GuCi/
It produces a plot of the gc content of the whole genome. It works on small genomes, like bacteria.
Hope you like!

At least state that you're the developer of this kinda-related tool that absolutely DOES NOT generate the kind of plot shown in the question.


text:
What is happening in the Zuker-Algorithm
I am trying to understand the Zuker Algorithm.
What is happening during initialisation? I know in Nussinov the diagonal is filled with 0s. Are the diagonals filled with infinite energy during Zuker´s initialisation?
Also, what is happening during recursion?
I know unlike Nussinov, Zuker pays attention to stacking energies and loop structures, not individual basepairs. But how? How does the algorithm calculate the energies for stacking and loop structures at position 1 if there is no info about its neighbors yet?
I understand there are 2 matrices - V and W. Is it correct that V is part of W? And is it correct that V keeps information about stack interactions and how much energy a stack of base pairs contributes - like keeping score on a multiplier?
Is there a code or example (I would prefer R), which demonstrates how the algorithm goes over a sequence? Like a demonstration. I thought of one myself, but I am not sure if its correct:
sequence <- c(1, 2, 3, 4, 5,6,7,8,9,10,11,12,13,14,15)

n <- length(sequence)

for (l in 2:n) {
   for (i in 1:(n - l + 1)) {
      j <- i + l - 1
      cat("i =", i, "j =", j, "\n")
   }
}

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.

I am using RNAFold (http://rna.tbi.univie.ac.at/cgi-bin/RNAWebSuite/RNAfold.cgi). I am suppose to write a shor text after my work and I am having diffuculties understanding certain espects of Zuker`s Algorithm. Nussinov is easy to understand and visualize on paper. And i was wondering how would implement Zuker in R or Python for example or how one would fill a dotplot using Zuker.

May i ask.. Whats the graph theory code in Python. I found some Python code on Github which is suppose to simulate Zuker, but if its not the correct Zuker Algorithm, then I would like to know.

Graph theory code for RNA structure is a separate question. This question is about Zuker's code - I don't know it - I thought Zuker used graph theory (which I know) but used matrix algebra instead. It's easy to get hold of the code just download the src for MFold and especially RNAfold. RNAfold is in C.


text:
Running bwa mem and samtools in one script on a linux cluster
I am trying to write a script that will run this line from https://github.com/kmceres/Mbovis_pangenome/blob/main/Parameter_file.md:
bwa mem -t 14 /path/to/contig/contigs.fasta -r contig_1P.fq.gz contig_2P.fq.gz |
   samtools view -Sb |
      samtools sort -@14 -o /path/to/contig/contig.sorted.bam

All of the files are in directories under their sequence number within a bwa_output (from when I ran bwa index) directory (e.g. /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992)
Here is my script:
#!/bin/bash

#SBATCH --job-name=bwa_mem
#SBATCH --account=lilianasalvador
#SBATCH --partition=standard
#SBATCH --nodes=1
#SBATCH --ntasks=1
#SBATCH --mem-per-cpu=16gb
#SBATCH --time=99:00:00

module load bwa/0.7.17
module load samtools/1.19.2

for file in *_contigs.fasta; do
    sample_name=${file%_contigs.fasta}
mkdir -p /groups/lilianasalvador/noah_100/bwa/bwa_mem_output/$sample_name

bwa mem -t 14 *_contigs.fasta /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR* -r *_1.fastq.00.0_0.cor.fastq.gz *_2.fastq.00.0_0.cor.fastq.gz | samtools view -Sb | samtools sort -@14 -o /groups/lilianasalvador/noah_100/bwa/bwa_mem_output

done

I'm getting this error:
Usage: bwa mem [options] <idxbase> <in1.fq> [in2.fq]
...
[main_samview] fail to read the header from "-".
[W::hts_set_opt] Cannot change block size for this format
samtools sort: failed to read header from "-"

I am relatively new to all of this and would appreciate some help editing my script
update 7/1: I tried running this
module load bwa/0.7.17
module load samtools/1.19.2

outBamsDir=/groups/lilianasalvador/noah_100/bwa/bwa_mem_output/

for file in *_1.fastq.00.0_0.cor.fastq.gz; do
    sample_name=${file%_1*}
file2=${file/_1/_2}
    mkdir -p "$outBamsDir"/"$sample_name"
    bwa mem \
      -t 14 \
/groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992/SRR8599992_contigs.fasta \
      "$file" "$file2" |
         samtools sort -@14 \
          -o "$outBamsDir/$sample_name/$sample_name".bam
done

and I got this error
/var/spool/slurm/d/job10075741/slurm_script: line 40: -t: command not found
/var/spool/slurm/d/job10075741/slurm_script: line 42: /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992/SRR8599992_contigs.fasta: Permission denied
/var/spool/slurm/d/job10075741/slurm_script: line 44: SRR8599992_1.fastq.00.0_0.cor.fastq.gz: command not found
[E::hts_open_format] Failed to open file " " : No such file or directory
samtools sort: can't open " ": No such file or directory
/var/spool/slurm/d/job10075741/slurm_script: line 48: -o: command not found

I also tried this:
module load bwa/0.7.17
module load samtools/1.19.2

for file in *_contigs.fasta; do
    sample_name=${file%_contigs.fasta}
mkdir -p /groups/lilianasalvador/noah_100/bwa/bwa_mem_output/$sample_name

    bwa mem -t 14 /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992/SRR8599992_contigs.fasta -r SRR8599992_1.fastq.00.0_0.cor.fastq.gz SRR8599992_2.fastq.00.1_0.cor.fastq.gz | samtools view -Sb | samtools sort -@14 -o /groups/lilianasalvador/noah_100/bwa/bwa_mem_output/SRR8599992

done

and got this error
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 648167 sequences (95550856 bp)...
[M::mem_process_seqs] Processed 648167 reads in 28.660 CPU sec, 28.709 real sec
[main] Version: 0.7.17-r1188
[main] CMD: bwa mem -t 14 -r SRR8599992_1.fastq.00.0_0.cor.fastq.gz /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992/SRR8599992_contigs.fasta SRR8599992_2.fastq.00.1_0.cor.fastq.gz
[main] Real time: 42.903 sec; CPU: 30.392 sec
[bam_sort_core] merging from 0 files and 14 in-memory blocks...
[E::hts_open_format] Failed to open file "/groups/lilianasalvador/noah_100/bwa/bwa_mem_output" : Is a directory
samtools sort: failed to create "/groups/lilianasalvador/noah_100/bwa/bwa_mem_output": Is a directory

This is what worked:
module load bwa/0.7.17
module load samtools/1.19.2

bwa mem -t 14 /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR8599992/SRR8599992_contigs.fasta -r SRR8599992_1.fastq.00.0_0.cor.fastq.gz SRR8599992_2.fastq.00.1_0.cor.fastq.gz |
   samtools view -Sb |
    samtools sort -@14 -o /groups/lilianasalvador/noah_100/bwa/bwa_mem_output/SRR8599992/SRR8599992.contig.sorted.bam

What exactly are you trying to do? Do you have one set of fastq files and you want to align that against multiple genomes (contigs) or do you have many sets of fastq files, the output of multiple sequencing runs (or the same run but multiple sequenced samples) and you want to align each set of fastqs against the contigs?

Please [edit] your question and explain what your files are, and then also show us the output of ls *_contigs.fasta, ls /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR*, ls *_1.fastq.00.0_0.cor.fastq.gz and ls  *_2.fastq.00.0_0.cor.fastq.gz. I tried to answer but without knowing what the actual data and end objective are, I can only make guesses.

Please read https://bioinformatics.stackexchange.com/help/formatting and format your code as code.

I see two probable sources of error. The general format of the bwa command you are trying to run is
bwa mem genome.fa read1.fastq.gz read2.fastq.gz

However, you are running (with some options removed and paths truncated for brevity):
bwa mem *_contigs.fasta /groups/[...]/SRR* -r *_1.[...].fastq.gz *_2.[..].fastq.gz

Now, I don't know what the globs *_contigs.fasta, /groups/[...]/SRR*, *_1.[...].fastq.gz or *_2.[..].fastq.gz would expand to but I doubt the result of expanding so many globs will be just one genome file and two fastq files. In addition, you seem to be already running this once for every file ending in _contigs.fasta since you have it in the loop:
for file in *_contigs.fasta; do ...; done

Assuming you actually have one file wth all of your contigs and this file has been indexed and that's what you want to use as a genome, then you would want
bwa mem /path/to/contig/contigs.fasta read1.fastq.gz read2.fastq.gz

I don't quite understand what you have from your question though so I can't be more specific. Next, presumably you have various sequencing run outputs, so multiple pairs of fastq files. Those are the ones your want to iterate over, not the contigs. So you would want a loop like this (I removed the unnecessary samtools view step, since samtools sort can handle the sam to bam conversion, and the also unnecessary -r option; if you know you need the -r, you can put it back):
outBamsDir=/groups/lilianasalvador/noah_100/bwa/bwa_mem_output/
for file in *_1.fastq.00.0_0.cor.fastq.gz; do
    sample_name=${file%_1*}
    file2=${file/_1/_2}
    mkdir -p "$outBamsDir"/"$sample_name"
    bwa mem \
      -t 14 \
      /path/to/contig/contigs.fasta \
      "$file" "$file2" |
         samtools sort -@14 \
          -o "$outBamsDir/$sample_name/$sample_name".bam
done

Or, if you really have multiple contig files and want to treat each of them as the reference genome to align against, so presumably you only have 2 fastq files, then the loop would be:
outBamsDir=/groups/lilianasalvador/noah_100/bwa/bwa_mem_output/
for file in *_contigs.fasta; do
  sample_name=${file%_contigs.fasta}
  mkdir -p "$outBamsDir"/"$sample_name"
  bwa mem -t 14 /whatever/you/want/to/use/as/genome.fa \
   "$file" |
   samtools sort -@14 -o "$outBamsDir/$sample_name".bam
done

That, of course, assumes single-end sequencing with one file per sample.

Good points. AFAIK neither bwa  nor bwa-mem2 can take more than 2 FASTQ files, guess the paired or non-paired mode and align them. Besides, using a single computing node to align multiple FASTQ pairs sequentially is not optimal save for a single node "cluster".

I have a contigs.fasta file for each sequence as well as two .cor.fastq.gz files (reverse and forward) for each sequence. In /groups/lilianasalvador/noah_100/bwa/bwa_output I have 100 files - each corresponding to a sequence (i.e. SRR8599992, SRR8599993...) and within each of those is the contigs.fasta and the .cor.fastq.gz files for the sequence. That's why I was saying _contigs.fasta /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR to try to explain that the contigs.fasta files are within the corresponding directory that starts with SRR. Idk how to correctly write the path

Do I have too many files? I was told I should have a contigs.fasta file for every sequence, which is why they're separated in directories

You might or you.might not, but we cannot know since you never showed us the files or even explained what you are really trying to do. However if you found a solution, please post that as an answer and not as part of the question.

I note that you make a mistake in your command lines,it seems that you want to map reads to the ref , convert the sam to the bam file and sort the bam file in one command, actually, you need to get the unsorted bam file first and sort the bam file with samtools singly. Because samtools sort will sort the whole file,but if you use the pipe as input of the 'samtools sort', the result is not a whole file from the samtools view -Sb,it just is a part of the resulty,so the samtools sort can't sort it , and the result may be stored in the memory , when many parts of results into the memory,it will be killed as OOM. You can try to map and convert sam to bam first and do the sorting singly.
I find that it contains /groups/lilianasalvador/noah_100/bwa/bwa_output/SRR between the reference sequence and the -r , it seems don't need this parameter , you can confirm the manual of the bwa mem

This isn't true. Samtools sort is perfectly capable of sorting from a pipe. Try it, just run samtools view -h foo.bam | samtools sort -o bar.bam. And it can do this with humongous files, it doesn't store in memory, it creates temporary files while sorting precisely to avoid the problem you mention.

Sorry about that,actually i have tried many times in my machine, but that's not hpc,it just has 64GB memory, and I indeed got that problem,when i use the bwa mem | samtools view | samtools sort my computer will be down, i will can't connect it by ssh, and after long time ,it will recover,but that is uncontrollable ,and i hava find out some oom events during the analysis when i can connect it again. I have resolved it by the way i hava said.But there is an other question in his commandline which has been mentioned in the second paragraph in my answer last post.It seems that you mention it too.

That can be caused by various things, but I can assure you, I have been running bwa mem | samtools sort in a professional, production pipeline for almost a decade now. And also did so on my laptop, before posting, just to double check. I also confirmed in the [samtools manual(https://www.htslib.org/doc/samtools-sort.html) that it creates temp files when memory is not enough. As for the last paragraph, I'm sorry but I don't understand it at all. Yes, the -r isn't needed, but it isn't causing the problem either.

Thanks to reply, actually i have noted the directory which he should not use between the reference and the reads file too. But the words may be unclear.  On the other hand i major in computer science ,when i get the  problem i have mentioned and take account for the behaviour of my comuputer and the syslog, i think i could be the reason what i have said , and by split the command line i resolve it. So it makes me believe the reason more.  The version of the samtools and htslib is 1.20. I will try to do reproduction experiments to located the ture reason. Thank you for reminding me.


text:
Extract gene sequences from isolates
I have 500 M.Bovis isolates and I want to extract the rpoB and katG genes from these sequences. I took the nucleotide sequences from NCBI and saved them as a fasta file, and then tried to blast for one isolate using:
blastn -db databases/reference \
  -query sequences.fasta \
  -evalue 1e-3 \
  -word_size 11 \
  -outfmt 0 > sequences.reference.

But my output contains other quality data apart from the sequences. Is there anyway to extract just the sequences without other information?

Thats easy simply use minimap
minimap2 --heap-sort=yes --splice-flank=no mygenomes.fastq.gz fastaDB.fas

where

mygenomes.fastq.gz

is the 500 M. bovis genomes and ...

fastaDB.fas

are the genes you are interrogating the database with
The default parameters are set for the problem you are facing, i.e. to extract target genes within a species. In addition, you can just load fastaDB.fas with all the genes you are searching and it will loop through all of them on a single run.
It does have complex parameterisation, but thats usually to lower the criteria for a hit between the query sequence and genomes database.
minimap can ingest directly from a zip file (big advantage), has some parallelisation abilities (3 CPUs used) and its very fast (faster than Blast).

Appreciate the help


text:
How to download multiple UniRef clusters with all their members?
We are trying to download multiple UniRef clusters in one batch on the website. However, it seems the members are not downloaded... Could you guide me a way to the solution? Many thanks!
For example, see this link. We would like to download all these clusters including their members. But when selected these clusters, only the representative sequences (one sequence per cluster) were downloaded.

Please [edit] your question and give us some context. What website? How are you doing this? Which clusters? What do you want to download exactly, the protein sequence?

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.

Please see the link below. We would like to download all these clusters including their members. But when selected these clusters, only the representative sequences (one sequence per cluster) were downloaded...
https://www.uniprot.org/uniref?query=membrane+protein&facets=identity%3A0.5

Can you explain why you want this? Are you sure UniRef is the right tool for what you are trying to do? How will the sequences be used? The whole point of UniRef is to give one representative sequence so it seems strange that you want to do something else. Please [edit] your question and explain what the final objective here is. Perhaps you are trying to get orthologous sequence groups?

The UniProt page apparently only lets you download fasta sequences for the representative protein in each UniRef group:

So you can't get all the fasta sequences from there. However, once you get the list of UniProt IDs, you can get the sequences from those. To get the IDs, first select TSV output:

Then remove everything and add "cluster members" from the "Miscellaneous" section:

This will give you a file that looks like this:
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz
Cluster ID  Cluster members
UniRef50_A0A0F6EQC6 A0A0F6EQC6; E3E6D7; A0A3G8RAV1; A0A2K1EFR7; A0A2S6NYI6; A0A7H0Y5M1; A0A2N9YZ66; A0A853Q3J3; A0A8I2C241; A0A0K2FC93
UniRef50_A0A0T6A5I5 A0A0T6A5I5; A0A2V7BZB0; A0A2V6XWF2
UniRef50_A0A2H4TXU1 A0A2H4TXU1
UniRef50_A0A6C2YI31 A0A6C2YI31
UniRef50_B1KXG0 B1KXG0; UPI001F3F1379

It's a tab separated file, and the second field is a ;-separated list of cluster members. Get those into a new file, one per line:
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz |
    tail -n +2 | cut -f2 | sed 's/; /\n/g'
A0A0F6EQC6
E3E6D7
A0A3G8RAV1
A0A2K1EFR7
A0A2S6NYI6
A0A7H0Y5M1
A0A2N9YZ66
A0A853Q3J3
A0A8I2C241
A0A0K2FC93
A0A0T6A5I5
A0A2V7BZB0
A0A2V6XWF2
A0A2H4TXU1
A0A6C2YI31
B1KXG0
UPI001F3F1379

Save to a file with
$ zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz |
    tail -n +2 | cut -f2 | sed 's/; /\n/g > newfile

It may be possible for the same protein to belong to multiple clusters, I don't know UniRef well enough to be sure if this is the case. So, to be on the safe side, make the IDs unique:
zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz | tail -n +2 | cut -f2 | sed 's/; /\n/g' | sort -u

That assumes GNU sort, the default on Linux. If yours doesn't have a -u option, use this instead:
zcat uniref_id_UniRef50_A0A0F6EQC6_OR_id_U_2024_03_01.tsv.gz | tail -n +2 | cut -f2 | sed 's/; /\n/g' | sort | uniq

Now, you can use these IDs and UniProt's ID mapping tool to get their sequences. Go to the mapping tool and paste in your IDs (this assumes you have no more than 100000 IDs. If you have more, do this in multiple steps.

Click on "Map N IDs" and that will take you to a page like this, where you can now click on "Download" to retrieve the sequences:

The result of all this is a fasta file like this:
$ zcat foo.fa.gz | head
>tr|A0A0F6EQC6|A0A0F6EQC6_PAEPO Membrane protease regulatory membrane protein OS=Paenibacillus polymyxa OX=1406 GN=yuaF PE=4 SV=1
METLFWWLLAGGAIFTVVSVLIGDVLGSWLDGALELPGLDWFRPVVLLGAMTAFGGAGVL
LTKYTGLTMGWVVVLALGIALVIGVLVFFAFIKPMANSEVSSGFSMRELTGRIGEVTIPV
PEVGYGEVMIRLGAGNTIHTASSFDHKPLAAGTRIVVVEVAEGVVRVAYLDT
>tr|A0A0K2FC93|A0A0K2FC93_9BACL Membrane protease regulatory membrane protein OS=Paenibacillus peoriae OX=59893 GN=ABE82_18870 PE=4 SV=1
METIFWWLLAGGAIFTVVSVLIGDVLGGWLDGALELPGLDWFRPVVLLGAMTAFGGAGVL
LTKYTGLSMGWVVVLALGIAFVIGILVFFAFIKPMANSEVSSGFSMRELTGRIGEVTIPV
PEVGYGEVMIRLGAGNTIHTASSFDHKPLAAGTRIVVVEVAEGVVRVAYLDT
>tr|A0A0T6A5I5|A0A0T6A5I5_9BACT Membrane protein, putative membrane protein OS=candidate division NC10 bacterium CSP1-5 OX=1640516 GN=XU15_C0006G0101 PE=4 SV=1
MDSFVQLVVGTVVLRPYVVAFLLLYLFAAVADMGWRRTLLFTVVAWSVAFTAESLSIRFG

However, I would urge you to think if this is really the best way for what you want to do. If you are trying to find orthologous groups, there are better ways.


text:
Databases for mtDNA with location/age
I'm doing some experiments on mtDNA, trying to infer samples' location/age based on the similarity among sequences. I was able to find good samples (with geographical coordinates and age range) for my model on https://amtdb.org/sample/, but other databases (e.g., https://www.ncbi.nlm.nih.gov/, https://www.mitomap.org/) don't look to have the desired data available.
Would you know some database, website, project, etc. that has mtDNA sequences available along with location and/or age?
Thank you.

NCBI is accessible via their API, with Biopython Bio.Seq method SeqIO and SeqRecord using Genbank format results in objects with very rich features on each sample. Thus it should be possible to parse.

What species - are these human mtDNA samples? What do you mean by age? Given the low mutation rate of mitochondria, it's not going to change much over the course of a single person's life.


text:
How to download the same 18S rRNA gene region for multiple species?
I need 18S rRNA gene sequences for a wide variety of eukaryotic species.
I know that Ensembl has this gene sequence, but I don't understand what is the correct protocol to download the sequences. I already know  the basics of Ensembl's Rest API, but there are still a few things that I don't understand.
I only know the gene name (18S rRNA) and the list of taxa I want to select it from. But I do not know the region of the gene.
For example: I tried looking for '18srrna' on the Ensembl browser and I noticed that for the Fruit fly there are 2 sequences corresponding to this gene (Ensembl IDs: FBgn0085813 and FBgn0267521).

How do I know which sequence I should use to make an alignment with other species?
How can I batch download the same region for all the species of interest?
Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

I don't have a solution for question (1) but for question (2) I am planning on requesting from the API the output of the endpoint GET sequence/region/:species/:region. But then again the problem is that I would need the position of the region for each species, and it sounds unfeasible to do for more than a dozen of species.

Hi @Mirk. I'll answer over the weekend if no-one else does. You wrote, "eucariotic" do you mean eukaryotes?

Yes, sorry I did a spelling mistake but now I corrected it. Thanks for pointing it out.

I can suggest two ways to batch-download aligned sequences for the species of interest. Both ways are integral to the Ensembl database and toolsets.
But before that, I suggest FBgn0085802 (which is rRNA) instead of FBgn0085813 and FBgn0267521 (which are pseudogenes). You may check FlyBase database.

Metazoa Ensembl is preferable Since you are interested in Fruit fly 18S rRNA homologs. You query Ensembl ID in the Ensembl search box (at www.metazoa.ensembl.org) and click the popping link. 18SrRNA:CR41548 comes for FBgn0085802. On the left-hand side, Gene-based displays, select Sequence. under the Marked-up sequence, pick "BLAST this sequence".

Opening the BLAST/BLAT search automatically sets a DNA alignment query. You may search against the species of interest and add "18S" as Description. Later, run without changing any other parameter. It may take several minutes. Then, view results for hits found. You can investigate Sequence/s with the best alignment and closest genomic locations, then export data in FASTA or other formats.
While this assures you of the alignment and lets you investigate in detail the homologs, you may not batch-download particularly in this way.

As mentioned, BioMart Restful API is a way to batch-download sequences that you want. The website does only allow one search at a moment but editing the Perl script, which is generated according to your query, you may get sequences for multiple datasets. For setup, you should follow BioMart Perl API guideline cautiosly.

For your example query, set proper datasets, filters and attributes. Again for fruit fly, you may refer to:
https://www.metazoa.ensembl.org/info/data/biomart/biomart_perl_api.html
My example query seems like below in popping Perl script:
$query->setDataset("dyakuba_eg_gene");
$query->addFilter("with_dmelanogaster_eg_homolog", ["Only"]);
$query->addFilter("biotype", ["rRNA"]);
$query->addFilter("transcript_biotype", ["rRNA"]);
$query->addAttribute("ensembl_gene_id");
$query->addAttribute("ensembl_transcript_id");
$query->addAttribute("gene_biotype");
$query->addAttribute("description");
$query->addAttribute("transcript_biotype");

How do I know which sequence I should use to make an alignment with other species?

Use both. For fruit flies there should thousands of 18S sequences.

How can I batch download the same region for all the species of interest?

Using Blast (below)

Do you know any protocol/guide that goes into detail of how to download and select sequences to make phylogenies?

Sure - Blast ... align via muscle -super5 and build via IQ-TREE
EMBL I don't know - using REST API via a CLI is possible. The problem is that you could potentially 1000s of 18S genes but a simple search term can ignore everything that isn't part of a genome explicitly labelled "18S". Thats where a Blast style approach is useful.
However, for NCBI Blast they intentionally block access to CLI with >100 hits requested. This is a complete pain and really, really limits using Blast to extract out a gene remotely using a remote request script.
Basically the way to do it is to ditch the API and run this via the Blast NCBI web interface, or download the entire eukaryotic database and perform local Blast, MMseqs2 or Minimap2.
It's a bit bizarre I know, NCBI just doesn't like bots running its Blast service, but it gives more permission to use its web-service. The central advantage is that Blast will nicely remove the 18S gene from any given genome which can be downloaded as a single locus, without needing to download loads of genome stuff thats not needed. NCBI still does impose a limit on how much can be blasted via the web which sadly now is only 5000 (I'm sure it was much more like 50000) and its really not cool.
What you can do is use the Blast search "exclude/include" parameters to target the taxa of interest within the 5000 limit and using a tree of relevant eukaryotic evolution navigate the data sets. I assume you are looking for around 100000 sequences. Basically, once you've blasted you have the option of downloading just the region blasted and nothing else ... or the whole genome (avoid this option)... you will then have 5000 nicely 18S genes excised in silico
You can also set a lot other blast parameters, in particular only considering full length hits.


text:
Problem while mapping reads to mtDNA (SortSam)
I am trying to map MiSeq reads to a reference genome and extract mutations using MToolBox, which implements gsnap, GATK, Picard, and other tools. When running the tool with example data, there were no errors, so I believe the tool is installed correctly. However, when I used the data sequenced with MiSeq for this run, I encountered the following error. Specifically, the OUTPUT file (OUT.sam.bam) from SORTING OUT.sam FILES WITH PICARDTOOLS was not generated. How should I address this issue? Could it be a problem with the sequence data?
Could someone please help me?
##### EXECUTING READ MAPPING WITH MAPEXOME...

mapExome for sample pc-mtGXL, files found: pc-mtGXL.R1.fastq.gz pc-mtGXL.R2.fastq.gz
Mapping onto mtDNA... /lustre7/home/User/MToolBox-1.2.1/MToolBox/bin/gmap/bin/gsnap -D /lustre7/home/User/MToolBox-1.2.1/MToolBox/gmapdb/ --gunzip -d chrM -A sam --nofails --pairmax-dna=500 --query-unk-mismatch=1
--read-group-id=sample --read-group-name=sample --read-group-library=sample --read-group-platform=sample -n 1 -Q -O -t 8 pc-mtGXL.R1.fastq.gz pc-mtGXL.R2.fastq.gz > /home/User/MTDNA-prospective/pc_GXL/out/OUT_pc-mtGXL/outmt.sam 2> /home/User/MTDNA-prospective/pc_GXL/out/OUT_pc-mtGXL/logmt.txt Extracting FASTQ from SAM... Mapping onto complete human genome...single reads Mapping onto complete human genome...pair reads Reading Results... Filtering reads... Outfile saved on /home/User/MTDNA-prospective/pc_GXL/out/OUT_pc-mtGXL/OUT.sam. Done.

SAM files post-processing...

##### SORTING OUT.sam FILES WITH PICARDTOOLS...

[Thu May 16 12:23:08 JST 2024] net.sf.picard.sam.SortSam INPUT=OUT.sam OUTPUT=OUT.sam.bam SORT_ORDER=coordinate TMP_DIR=[/home/User/MTDNA-prospective/pc_GXL/out/OUT_pc-mtGXL/tmp]    VERBOSITY=INFO QUIET=false VALIDATION_STRINGENCY=STRICT COMPRESSION_LEVEL=5 MAX_RECORDS_IN_RAM=500000 CREATE_INDEX=false CREATE_MD5_FILE=false [Thu May 16 12:23:08 JST 2024] Executing as User@at138 on Linux 5.15.0-87-generic amd64; OpenJDK 64-Bit Server VM
1.8.0_392-b08; Picard version: 1.98(1547) [Thu May 16 12:23:08 JST 2024] net.sf.picard.sam.SortSam done. Elapsed time: 0.01 minutes. Runtime.totalMemory()=2058354688 To get help, see http://picard.sourceforge.net/index.shtml#GettingHelp Exception in thread "main" net.sf.samtools.SAMFormatException: SAM validation error: ERROR: Read name M02699:57:000000000-LG9CK:1:1119:16584:12259, CIGAR M operator maps off end of reference     at net.sf.samtools.SAMUtils.processValidationErrors(SAMUtils.java:448)  at net.sf.samtools.SAMRecord.getCigar(SAMRecord.java:606)   at net.sf.samtools.SAMRecord.getCigarLength(SAMRecord.java:617)     at net.sf.samtools.SAMRecord.isValid(SAMRecord.java:1599)   at net.sf.samtools.SAMLineParser.parseLine(SAMLineParser.java:328)  at net.sf.samtools.SAMTextReader$RecordIterator.parseLine(SAMTextReader.java:237)   at net.sf.samtools.SAMTextReader$RecordIterator.next(SAMTextReader.java:225)    at net.sf.samtools.SAMTextReader$RecordIterator.next(SAMTextReader.java:201)    at net.sf.samtools.SAMFileReader$AssertableIterator.next(SAMFileReader.java:672)    at net.sf.samtools.SAMFileReader$AssertableIterator.next(SAMFileReader.java:650)    at net.sf.picard.sam.SortSam.doWork(SortSam.java:68)    at net.sf.picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:177)   at net.sf.picard.cmdline.CommandLineProgram.instanceMainWithExit(CommandLineProgram.java:119)   at net.sf.picard.sam.SortSam.main(SortSam.java:57) Success.

samtools index: "OUT.sam.bam" is in a format that cannot be usefully indexed

##### REALIGNING KNOWN INDELS WITH GATK INDELREALIGNER...

Realigning known indels for file OUT_pc-mtGXL/OUT.sam.bam using /home/User/MToolBox-1.2.1/MToolBox/data/MITOMAP_HMTDB_known_indels.chrM as reference... INFO  12:23:11,978 HelpFormatter -
------------------------------------------------------------------------------------  INFO  12:23:11,984 HelpFormatter - The Genome Analysis Toolkit (GATK) v3.8-1-0-gf15c1c3ef, Compiled 2018/02/19 05:43:50  INFO  12:23:11,984 HelpFormatter - Copyright (c) 2010-2016 The Broad Institute  INFO  12:23:11,984 HelpFormatter - For support and documentation go to https://software.broadinstitute.org/gatk  INFO  12:23:11,984 HelpFormatter - [Thu May 16 12:23:11 JST 2024] Executing on Linux
5.15.0-87-generic amd64  INFO  12:23:11,984 HelpFormatter - OpenJDK 64-Bit Server VM 1.8.0_392-b08  INFO  12:23:11,986 HelpFormatter - Program Args: -U ALLOW_N_CIGAR_READS -T IndelRealigner -R /home/User/MToolBox-1.2.1/MToolBox//data/chrM.fa -I OUT.sam.bam -o OUT.realigned.bam -targetIntervals /home/User/MToolBox-1.2.1/MToolBox//data/intervals_file_chrM.list
-known /home/User/MToolBox-1.2.1/MToolBox//data/MITOMAP_HMTDB_known_indels_chrM.vcf
-compress 0  INFO  12:23:11,993 HelpFormatter - Executing as User@at138 on Linux 5.15.0-87-generic amd64; OpenJDK 64-Bit Server VM
1.8.0_392-b08.  INFO  12:23:11,994 HelpFormatter - Date/Time: 2024/05/16 12:23:11  INFO  12:23:11,994 HelpFormatter -
------------------------------------------------------------------------------------  INFO  12:23:11,994 HelpFormatter -
------------------------------------------------------------------------------------  INFO  12:23:12,139 NativeLibraryLoader - Loading libgkl_compression.so from jar:file:/lustre7/home/User/MToolBox-1.2.1/MToolBox/ext_tools/GenomeAnalysisTK.jar!/com/intel/gkl/native/libgkl_compression.so INFO  12:23:12,711 GenomeAnalysisEngine - Deflater: IntelDeflater  INFO  12:23:12,711 GenomeAnalysisEngine - Inflater: IntelInflater  INFO  12:23:12,711 GenomeAnalysisEngine - Strictness is SILENT  INFO  12:23:12,962 GenomeAnalysisEngine - Downsampling Settings: No downsampling  INFO  12:23:12,967 SAMDataSource$SAMReaders - Initializing SAMRecords in serial  INFO  12:23:12,985 SAMDataSource$SAMReaders - Done initializing BAM readers: total time
0.02
##### ERROR ------------------------------------------------------------------------------------------
##### ERROR A USER ERROR has occurred (version 3.8-1-0-gf15c1c3ef):
##### ERROR
##### ERROR This means that one or more arguments or inputs in your command are incorrect.
##### ERROR The error message below tells you what is the problem.
##### ERROR
##### ERROR If the problem is an invalid argument, please check the online documentation guide
##### ERROR (or rerun your command with --help) to view allowable command-line arguments for this tool.
##### ERROR
##### ERROR Visit our website and forum for extensive documentation and answers to
##### ERROR commonly asked questions >https://software.broadinstitute.org/gatk
##### ERROR
##### ERROR Please do NOT post this error to the GATK forum unless you have really tried to fix it yourself.
##### ERROR
##### ERROR MESSAGE: SAM/BAM/CRAM file OUT.sam.bam is malformed. >Please see https://software.broadinstitute.org/gatk/documentation/article?id=1317for more information. Error details: SAM file doesn't have any read groups defined in the header.  The GATK no longer supports SAM files without read groups
##### ERROR ------------------------------------------------------------------------------------------

The last process reported an error. Exit.

OUT.sam.bam is 0 KB (i.e. empty). It is human mtDNA.

Seems that MToolBox is quite frozen in time (last release (Sep 2019) with hg19 and recommended tools like GSNAP 2015-12-31.v7 and samtools 1.3. Are there any alternative packages/pipelines being maintained?

The first error seems to come from picard when sorting the sam file:

SAMFormatException: SAM validation error: ERROR: Read name M02699:57:000000000-LG9CK:1:1119:16584:12259, CIGAR M operator maps off end of reference

From picard FAQ:

Picard validation errors may be turned into warnings by passing the command line argument VALIDATION_STRINGENCY=LENIENT.
Picard validation messages may be suppressed completely with VALIDATION_STRINGENCY=SILENT.
Another option is to use the CleanSam tool to soft-clip these reads so they don't map off the end of the reference.


text:
BLAST Database error: No alias or index file found for nucleotide database
I have created a database and added environment variables, but an "error message" still appears. How can I solve it?
error message: BLAST Database error: No alias or index file found for nucleotide database [/tmp/ramdisk_prok/nt_prok] in search path [/home/Tina:/tmp/ramdisk_prok:]

Blast Database
$ cd /tmp/ramdisk_prok

    nt_prok.00.nhr  nt_prok.04.nnd  nt_prok.08.nog  nt_prok.13.nhr  nt_prok.17.nnd
    nt_prok.00.nin  nt_prok.04.nni  nt_prok.08.nsq  nt_prok.13.nin  nt_prok.17.nni
    nt_prok.00.nnd  nt_prok.04.nog  nt_prok.09.nhr  nt_prok.13.nnd  nt_prok.17.nog
    nt_prok.00.nni  nt_prok.04.nsq  nt_prok.09.nin  nt_prok.13.nni  nt_prok.17.nsq
    nt_prok.00.nog  nt_prok.05.nhr  nt_prok.09.nnd  nt_prok.13.nog  nt_prok.18.nhr
    nt_prok.00.nsq  nt_prok.05.nin  nt_prok.09.nni  nt_prok.13.nsq  nt_prok.18.nin
    nt_prok.01.nhr  nt_prok.05.nnd  nt_prok.09.nog  nt_prok.14.nhr  nt_prok.18.nnd
    nt_prok.01.nin  nt_prok.05.nni  nt_prok.09.nsq  nt_prok.14.nin  nt_prok.18.nni
    nt_prok.01.nnd  nt_prok.05.nog  nt_prok.10.nhr  nt_prok.14.nnd  nt_prok.18.nog
    nt_prok.01.nni  nt_prok.05.nsq  nt_prok.10.nin  nt_prok.14.nni  nt_prok.18.nsq
    nt_prok.01.nog  nt_prok.06.nhr  nt_prok.10.nnd  nt_prok.14.nog  nt_prok.19.nhr
    nt_prok.01.nsq  nt_prok.06.nin  nt_prok.10.nni  nt_prok.14.nsq  nt_prok.19.nin
    nt_prok.02.nhr  nt_prok.06.nnd  nt_prok.10.nog  nt_prok.15.nhr  nt_prok.19.nnd
    nt_prok.02.nin  nt_prok.06.nni  nt_prok.10.nsq  nt_prok.15.nin  nt_prok.19.nni
    nt_prok.02.nnd  nt_prok.06.nog  nt_prok.11.nhr  nt_prok.15.nnd  nt_prok.19.nog
    nt_prok.02.nni  nt_prok.06.nsq  nt_prok.11.nin  nt_prok.15.nni  nt_prok.19.nsq
    nt_prok.02.nog  nt_prok.07.nhr  nt_prok.11.nnd  nt_prok.15.nog  nt_prok.20.nhr
    nt_prok.02.nsq  nt_prok.07.nin  nt_prok.11.nni  nt_prok.15.nsq  nt_prok.20.nin
    nt_prok.03.nhr  nt_prok.07.nnd  nt_prok.11.nog  nt_prok.16.nhr  nt_prok.20.nnd
    nt_prok.03.nin  nt_prok.07.nni  nt_prok.11.nsq  nt_prok.16.nin  nt_prok.20.nni
    nt_prok.03.nnd  nt_prok.07.nog  nt_prok.12.nhr  nt_prok.16.nnd  nt_prok.20.nog
    nt_prok.03.nni  nt_prok.07.nsq  nt_prok.12.nin  nt_prok.16.nni  nt_prok.20.nsq
    nt_prok.03.nog  nt_prok.08.nhr  nt_prok.12.nnd  nt_prok.16.nog  taxdb.btd
    nt_prok.03.nsq  nt_prok.08.nin  nt_prok.12.nni  nt_prok.16.nsq  taxdb.bti
    nt_prok.04.nhr  nt_prok.08.nnd  nt_prok.12.nog  nt_prok.17.nhr  taxdb.tar
    nt_prok.04.nin  nt_prok.08.nni  nt_prok.12.nsq  nt_prok.17.nin  taxonomy4blast.sqlite3

$ less -S ~/.bashrc

    export PATH=$PATH:/staging/tsailab/Blast/ncbi-blast-2.15.0+/bin
export BLASTDB=/tmp/ramdisk
export PATH=$PATH:/tmp/ramdisk
    export PATH=$PATH:/tmp/ramdisk_prok
    export BLASTDB=/tmp/ramdisk_prok

$ blastdbcmd -info -db /tmp/ramdisk_prok/nt_prok

    BLAST Database error: No alias or index file found for nucleotide database [/tmp/ramdisk_prok/nt_prok] in search path [/home/tsailab:/tmp/ramdisk_prok:]

$ blastn -query /staging/Blast/Input/HBV/Sample1.seq.fasta -db /tmp/ramdisk_prok/nt_prok -task blastn -evalue 1e-6 -max_hsps 1 -max_target_seqs 500 -outfmt "6 qseqid qlen qstart qend sseqid slen sstart send length nident pident mismatch gaps gapopen qcovs evalue bitscore stitle sskingdoms staxids sscinames sblastnames qseq sseq" -num_threads 10 -out /staging/Blast/Output/HBV/Sample1.original.txt*

    BLAST Database error: No alias or index file found for nucleotide database [/tmp/ramdisk_prok/nt_prok] in search path [/home/Tina:/tmp/ramdisk_prok:]

First off, why are you doing this in /tmp? Anything in /tmp will be deleted when the machine reboots. Next, are your DB files in /tmp/ramdisk_prok or in /tmp/ramdisk_prok/nt_prok? You seem to show that they're in /tmp/ramdisk_prok but then you give blast -db /tmp/ramdisk_prok/nt_prok, which is it?

The easiest solution is to install miniconda and then install via mamba
mamba create --name myenvname blast
mamba activate blast
mamba install blast
mamba update blast
blastn -query /staging/Blast/Input/HBV/Sample1.seq.fasta -db /tmp/ramdisk_prok/nt_prok -task blastn -evalue 1e-6 -max_hsps 1 -max_target_seqs 500 -outfmt "6 qseqid qlen qstart qend sseqid slen sstart send length nident pident mismatch gaps gapopen qcovs evalue bitscore stitle sskingdoms staxids sscinames sblastnames qseq sseq" -num_threads 10 -out /staging/Blast/Output/HBV/Sample1.original.txt
mamba deactivate

Check
Do you have io permission for
/staging/Blast/Input/HBV/Sample1.seq.fasta
/staging/Blast/Output/HBV/Sample1.original.txt

To install mamba you'll need Anaconda because its within the Anaconda collection and simplifies everything on three counts, viz:

in answer to your question all configurations are automated
it is easy to automate all updates
it offers a virtual environment where code can't mess with your system
whenever you want to download other packages its the same approach

Thus all you have to do is download and install Anaconda or miniconda. From this point everything is done using conda, either:
Anaconda installs in two different locations depending on whether you install the GUI version, which makes maintenance and environments a bit easier. In the later case I think its /opt, otherwise the commandline version is ~/anadonda3.
Thereon,
conda config --add channels conda-forge
conda update -n base --all
conda install -n base mamba


text:
Cannot run short-ranged nonbonded interactions on a GPU because no GPU is detected
I use these commands to create a build directory and configure the build:
mkdir build
cd build
cmake .. -DGMX_GPU=CUDA \
         -DCMAKE_CUDA_COMPILER=/usr/local/cuda/bin/nvcc \
         -DCMAKE_INSTALL_PREFIX=/usr/local/gromacs \
         -DGMX_BUILD_OWN_FFTW=ON -DGMX_MPI=OFF

and when I run this command for running a simulation
gmx mdrun -deffnm nvt -nb gpu

it is giving me this error:
Fatal error:
Cannot run short-ranged nonbonded interactions on a GPU because no GPU is detected.

please help me if you know.

Please [edit] your question and give us (a lot) more context. First of all, what is your operating system? Next, what are you building? Where did you download it from?  How can we try to reproduce this? What is your hardware? Do you have an NVIDIA graphics card? Is it compatible with whatever you are using? Are those the only commands you ran? Did you skip make install and sourcing /usr/local/gromacs/bin/GMXRC as described in the documentation?


text:
Extracting RNA reads alignments from a mix of DNA and RNA alignments
I have a fastq file that somehow has a mix of DNA and RNA. It's the sequence data of a DNA sample, but I found one interesting gene that has much higher coverage in exons than in introns. I used STAR to re-perform the alignment, and saw the clear exon-intron-exon pattern in the reads in IGV. So I suspect that there is RNA contamination in this sample.
This is whole-genome sequencing data. We used bwa at the beginning, and then STAR to check the exon-intron pattern. The library was prepared according the common DNA library kit's manual, as this was prepared by a company.
In order to estimate whether it is indeed a whole-genome RNA contamination in this sample and how much contamination there is in this sample, I want to estimate how many genes have RNA sequences aligned. I thought to use a splice-aware mapping tool to do the mapping, and then extract records that have gaps in them by filtering the CIGAR score for N. But I then realized that this may not be accurate, because alignment is performed against the reference genome and not the genome of the individual who was sequenced. If I filter the CIGAR score for N, I would also include the reads that bear large deletions in their own genome compared to the reference genome, even if it's DNA sequences.
So now I can't think of a way to screen the genes which have RNA sequences aligned. Without extracting the RNA reads, I also can't track where the RNA contamination come from, whether it's from our own sample or from the sequencing facility.
Another method I could think of is to filter the genes by their very different coverages of reads. But I also don't have a clear idea on how to do it.
So could anyone help me a bit here how to extract RNA reads or how to screen genes that have RNA reads aligned?

Some IGV screenshots would certainly help. Isolating reads from RNA could be tough, any chance you can resequence the sample?

@KarthikNair We're trying to avoid resequencing, because we actually have a whole population that has similar problems.

Couldn't you simply map against the transcriptome and then filter for alignments that are not split reads and of high quality? Essentially, if there is notable RNA contamination (hard to believe since RNA is fragile in the presence of RNases which you have abundantly when you break open a cell using a non RNA-tailored buffer) -- then you should have lots of exon-exon junctions covered.


text:
How can I visualize WDL workflows?
I would like to graphically visualize my workflow. I tried with Pipeline Builder but I didn't succeed as the http://pb.opensource.epam.com webpage doesn't work.
Thanks!

Putting your question title into a search engine returned https://github.com/miniwdl-ext/wdlviz, is that what you need? If not, please [edit] your question and describe your requirements. What do you want to give as input? What kind of output are you expecting? Can you give us a simple example to play with?


text:
calling SNPs on PacBio fltnc_mapped.bam reads
I am trying to iterate over PacBio fltnc_mapped.bam reads for SNPs. I am getting 92% filtering success with T-base SNP iteration and 74% filtering success with C-base SNP at the same position. All the reads  have MAPQ 60. Is there a way of stringently calling the SNPs or this is as good as it gets? My script:
# Load module
module load samtools

# Extract reads at SNP position
samtools view -h fltnc_mapped.bam chr1:138039992-138039992 > temp_reads.sam

# Separate the header and the reads for easier manipulation
grep "^@" temp_reads2.sam > header2.sam
grep -v "^@" temp_reads2.sam > reads2.sam

# Filter reads based on the SNPs
awk -v pos=138039992 -v ref_allele=T -v alt_allele=C '
  BEGIN { FS="\t"; OFS="\t" }
  /^@/ { print $0; next }  # Print header lines
  {
  # Calculate position in read corresponding to SNP
  for (i = 1; i <= length($10); i++) {
    read_pos = $4 + i - 1
    if (read_pos == pos) {
    base = substr($10, i, 1)
    if (base == ref_allele) {
    print $0 > "Ptprcb_reads2.sam"
    } else if (base == alt_allele) {
      print $0 > "Ptprca_reads2.sam"
    }
   }
  }
 }' temp_reads2.sam

# Add the header back to the filtered reads SAM files
cat header2.sam Ptprca_reads2.sam > Ptprca_reads_with_header2.sam
cat header2.sam Ptprcb_reads2.sam > Ptprcb_reads_with_header2.sam

# Convert SAM files to BAM files
samtools view -bS Ptprca_reads_with_header2.sam > Ptprca_reads2.bam
samtools view -bS Ptprcb_reads_with_header2.sam > Ptprcb_reads2.bam

# Clean up intermediate files
rm reads2.sam header2.sam Ptprca_reads2.sam Ptprcb_reads2.sam
  Ptprca_reads_with_header2.sam Ptprcb_reads_with_header2.sam

There are definitely better options than a simple mapped consensus, especially if you're using a low[er] accuracy / high-speed mapper like minimap2.
Have a look at Clair3; that seems to be fairly popular at the moment for long reads:
https://github.com/HKU-BAL/Clair3


text:
Keggscape Network is not showing nodes properly
I need to know how to fix this issue. No matter what pathway I use, the nodes are not displayed as shown on the keggscape website. Instead, they show up like this:

For reference, I used this pathway: https://www.genome.jp/pathway/hsa04152#

Please [edit] your question and highlight the specific issue. Is it that the node names aren't within the node boxes? Or is it something else?

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.

please refer to the link to see how the network is supposed to be like for the best understanding of the issue thank you

We are aware of the issue and are working on it. For some reason, the Keggscape styles are not getting applied but it is now fixed.  Please install the latest version.

Can you kindly described the versions where the bug was discovered and the version where it was fixed? The reason is that each post exists long after posting as a reference to others


text:
Finding gene duplications in bacteria
I would like to know if there is an easy tool to find how many copies of a certain gene are present in all available bacterial genomes.
For example, which bacterial genomes contain one or more copies of Superoxide dismutase.

Simply use RNAseq tools such as EdgeR, each gene must be a separate entity, e.g. separate fasta entry for this approach. This is one approach, there are lots. Probably a good idea to provide more details.

What is the dataset you are working from? Are you working from finished, annotated genomes, raw sequencing data, de novo assemblies, ...?
The simplest answer would be to go to a genome database such as KEGG and look up the relevant genes. MicrobesOnline and COGs are other resource I've used, but there are likely newer ones.
If you simply query the gene in such a database, you should get a list of entries and the relevant genomes that they're in, and you then get the copy number of the gene in the genome for free. Here is such a page on MicrobesOnline.
The one caution that I have is to really understand the orthology of the gene you're interested in.

Prokka is an easy-to-use pipeline for the annotation of bacterial genomes. Here is the link to the GitHub page: https://github.com/tseemann/prokka
nic


text:
What do the different node colors mean in the Cytoscape ReactomeFIViz functional interaction network?
I'm a bit stuck trying to figure out how Cytoscape ReactomeFIViz decides how nodes are colored in my Cytoscape functional interaction network generated from my input RNAseq dataset (table).
Can you please clarify?

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.

Could you show us an image of what you are seeing, perhaps? That could help. Have you looked at any of the tool's documentation?

Thank you terdon! I am unable to find the share/upload image option. However, what I am seeing looks something like the image under the 'Cluster FI network' paragraph on reactome.org/tools/reactome-fiviz#Gene_Set.2FMutation_Analysis . I looked through this documentation & under 4.1 Gene Set/Mutation Analysis
it says
Cluster FI network: run a network clustering algorithm (spectral partition based network clustering by Newman 2006) on the displayed FI network. Nodes in different network modules

will be shown in different colors (different colors used only for first 15 modules based on sizes).
What does this mean exactly?
Does the same apply to my analysis which uses
4.3 Microarray Data Analysis (even though the data I am analyzing is RNAseq)
?

I believe the clustering is based on connectivity, not your fold change values, so it would matter exactly how you got your network.