text:
arm64 installation of AMRFinderPlus
I was trying to install NCBI AMRFinderPLUS on my Nvidia Jetson using standard Bioconda based installtion. However, all the time it is giving error that the package is not available and I thought I am making some typo or channel is not proper.
(amrfinder) jay@jetmimd-desktop:~$ conda install -y -c conda-forge -c bioconda ncbi-amrfinderplus
Channels:
 - conda-forge
 - bioconda
 - defaults
Platform: linux-aarch64
Collecting package metadata (repodata.json): done
Solving environment: failed

PackagesNotFoundError: The following packages are not available from current channels:

  - ncbi-amrfinderplus

Current channels:

  - https://conda.anaconda.org/conda-forge
  - https://conda.anaconda.org/bioconda
  - defaults

To search for alternate channels that may provide the conda package you're
looking for, navigate to

    https://anaconda.org

and use the search bar at the top of the page.

Then I remembered about the different architecture. So, now here I am. Is there any way to install this on arm64 architecture?
Thank you for the help in advance!!!!

Try the following
conda install bioconda::ncbi-amrfinderplus

The above does work I tested it (below). There could be an issue with you version of conda (24.1.2 OSX).
Also has the package has been ported for GPU(?) ... I assume it is, otherwise it's obviously not going to work and conda will not install it
However, if it ported for GPU and still fails use mamba (I don't use mamba BTW), because it is very slow and that could be where the error is (unusually slow).
My connection initially was unsuccessful but resolved via

Retrying with flexible solve.

To demonstrate I got it working here's the reply ...
The following NEW packages will be INSTALLED:

  blast              bioconda/osx-64::blast-2.12.0-pl5262h78c34c6_0
  c-ares             pkgs/main/osx-64::c-ares-1.19.1-h6c40b1e_0
  curl               pkgs/main/osx-64::curl-7.88.1-hdb2fb19_2
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  hmmer              bioconda/osx-64::hmmer-3.3.2-hb19a6fa_3
  libcurl            pkgs/main/osx-64::libcurl-7.88.1-ha585b31_2
  libev              pkgs/main/osx-64::libev-4.33-h9ed2024_1
  libidn2            pkgs/main/osx-64::libidn2-2.3.4-h6c40b1e_0
  libnghttp2         pkgs/main/osx-64::libnghttp2-1.52.0-h1c88b7d_1
  libssh2            pkgs/main/osx-64::libssh2-1.10.0-hdb2fb19_2
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  ncbi-amrfinderplus bioconda/osx-64::ncbi-amrfinderplus-3.11.14-hfee4e81_0
  perl               pkgs/main/osx-64::perl-5.26.2-h4e221da_0
  perl-archive-tar   bioconda/osx-64::perl-archive-tar-2.32-pl526_0
  perl-carp          bioconda/osx-64::perl-carp-1.38-pl526_3
  perl-common-sense  bioconda/osx-64::perl-common-sense-3.74-pl526_2
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.087-pl526h6de7cb9_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.087-pl526h770b8ee_0
  perl-exporter      bioconda/osx-64::perl-exporter-5.72-pl526_1
  perl-exporter-tiny bioconda/osx-64::perl-exporter-tiny-1.002001-pl526_0
  perl-extutils-mak~ bioconda/osx-64::perl-extutils-makemaker-7.36-pl526_1
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.087-pl526h6de7cb9_0
  perl-io-zlib       bioconda/osx-64::perl-io-zlib-1.10-pl526_2
  perl-json          bioconda/osx-64::perl-json-4.02-pl526_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl526h04f5b5a_3
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-0.428-pl526_1
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.428-pl526_0
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl526h1de35cc_1
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.52-pl526h01d97ff_0
  perl-types-serial~ bioconda/osx-64::perl-types-serialiser-1.0-pl526_2
  perl-xsloader      bioconda/osx-64::perl-xsloader-0.24-pl526_0
  wget               pkgs/main/osx-64::wget-1.21.4-ha585b31_1

There is a blast gpu and I'd assume the other library/packages have gpu versions. There's a lot of Perl code there BTW, so it's an old package, before NVIDIA was big perhaps?

If there is a GPU compliant channel on conda (nanoporetch) - that is impressive: especially given the amount of Perl code. Porting for GPU from CPU isn't trivial.

Few versions behind

The key thing is whether the upgrades are fixing bugs or refactoring/adding functionality. If it's the latter then I would go  for it because the GPU would compensate for the slower runtime (obviously given you don't need the extra functionality).

Hey, Thank you very much! I had tried that yesterday  (among the many other things I could try/knew how to try) after posting here. But that also generated extact same error as posted above. I think(it may be basic and I am still not sure) the packages which work for OSX64 may not work for aarch64/arm64 and that is why I am not able to use this streamline way of installation. Although I also came across, separate channel on conda nanoporetch who have stated many packages for aarch64, just one or two versions behind. So, I am still in dilemna what to do!  Use the older versions or ...

infact I am not able to install the packages with noarch tag. Now, I am not sure weather this is a trivial information to know how to work around these packages or this is something which will take time in developing.

@prekij if there's no available package, you will probably have to compile it yourself for the architecture :/


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.

source: Stack Exchange

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

sort: votes


text:
How do I quantify a specific somatic variant?
I'm working with targeted Illumina sequencing data generated with DNA from diseased and healthy tissue (this is an age-related disease and is not cancer/neoplastic). My hypothesis is that the diseased tissue might contain specific somatic mutations in one particular gene (let's call it GENE), while the healthy tissue will not. The reasoning behind this: The same mutations cause autosomal-dominant versions of my disease of interest when present in the germline (my diseased tissue is from sporadic disease). My coverage varies, but is in the range of 1000x - 50'000x. My sequencing protocol employed UMIs.
I have run Mutect2 in default settings and with lower thresholds (--tumor-lod-to-emit) and haven't found any calls that correspond to these pathogenic mutations, so my results are negative. Nevertheless, I have reason to believe that my pathogenic mutations might be present at very low VAF (around 0.05-0.1 %), based on other experiments.
My question: How do I quantify the number of reads that correspond to a particular variant of interest? Let's say that base 2000 in GENE is usually an A, whereas the mutation A>T would be pathogenic. How do I count the number of Ts in position 2000?
As I'm somewhat new to bioinformatics, I'm also very open to any other ideas. Thanks!

Edit: @terdon's approach worked perfectly in cases where I used a VCF for --alleles that had been generated as an output of Mutect2, such as the following one:
VCF produced by Mutect:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Curiously, if I used a manually produced file such as the following, no output was produced, unless the variant of interest was heterozygous:
Manually produced VCF:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT
chr20   4699605 rs1799990   A   G   .   .   .

Output using the manual VCF, heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  heterogyzous_sample
chr20   4699605 .   A   G   .   .   AS_SB_TABLE=881,932|920,847;DP=3791;ECNT=1;MBQ=20,17;MFRL=248,243;MMQ=60,60;MPOS=21;POPAF=7.30;TLOD=2014.57 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1813,1767:0.469:3580:559,218:384,191:1720,1679:881,932,920,847

Output using the manual VCF, non-heterozygous sample:
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  novariant_sample

Manual inspection of this non-heterozygous sample revealed that it contained 1 read with the ALT variant, so this doesn't seem to depend on the variant being absent.

I have a little function that can report this by using samtools to look at the read pileup. Um. It isn't the most elegant piece of code, but it does sorta do the job:
variantSupport ()
{
  genome=$1
  bam=$2;
  range=$3;
  samtools mpileup -f "$genome" -r "$range" "$bam" | cut -f 5 |
    tr '[a-z]' '[A-Z]' | fold -w 1 | sort | uniq -c |
    awk 'BEGIN{tot=1}$2~/,/{ $2="." }
         {
           tot+=$1;
           c[$2]+=$1
         }
         $2~/[a-zA-Z]/{ varReads+=$1; }
         END{
           for(i in c){
             print i,c[i]
           }
           printf "%s of %s reads (%.1f%%), support the variant.\n",varReads,tot,varReads*100/tot
   }'
}

Add the function to your ~/.bashrc or just paste it directly into a terminal and you can then do:
$ variantSupport /path/to/genome/hg19.fa foo.bam chr13:28607916-28607916
[mpileup] 1 samples in 1 input files
A 915
* 1
. 93
915 of 1010 reads (90.6%), support the variant.

This is telling me that at position 28607916 of chr13, I have 1 read supporting a deletion (*, see https://www.htslib.org/doc/samtools-mpileup.html for an explanation of the symbols), 93 reads supporting the reference residue (T, in this case) and 915 supporting an A, so a T=> A variant.
However, a better approach is to use a variant caller and force it to call the variant you want. You can then investigate the support for it. With mutect2, you can use the --alleles option to force the calling and then use the -L option to limit calling to only the specific region. If you give the same vcf file for both options, that will return results for the variants you want and only those. For example, given this file foo.vcf.gz (which needs to be compressed using bgzip and indexed using tabix):
$ zcat 1.vcf.gz
##fileformat=VCFv4.3
##reference=/genomes/hg19/hg19.fa
##FORMAT=<ID=AB,Number=A,Type=Float,Description="Allelic Balance">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Read depth">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Net Genotype quality across all datasets, calculated from GQ scores of callsets supporting the consensus GT, using only one callset from each dataset">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##contig=<ID=chr13,assembly=hg19,length=115169878>
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  HG003
chr13   28607917    .   G   C   .   .   .   GT:AB:DP:GQ 1/1:0.6:395:318

You can then do:
gatk Mutect2 \
  -R /genomes/hg19/hg19.fa \
  --alleles 1.vcf.gz  \
  -L 1.vcf.gz \
  -I foo.bam \
  -O out.vcf

This will produce an out.vcf like this (I am excluding most headers for clarity):

##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  331970
chr13   28607917    .   G   C   .   .   DP=1107;ECNT=1;MBQ=30,0;MFRL=353,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-2.716e+00 GT:AD:AF:DP:F1R2:F2R1:SB    0/1:1043,0:9.616e-04:1043:518,0:496,0:989,54,0,0

I can see this is not a real variant since the value for AD, the number of reads supporting the reference and the variant is 1043,0, meaning that 1043 reads support the reference and none support the variant.
This way, you can get results for your variants of interest and see if there is any evidence supporting them.

Thanks! Your answer does what I needed. However, I have noticed that the Allele Depth numbers in the VCF (out.vcf, in your code) are lower than expected. For instance, I got REF=105 and ALT=104, whereas the numbers are  REF = 436 and ALT = 380 if I extract this information with bcftools and set -d to a high number. I need the total read numbers, so how can I edit your command to achieve this?

@Nereus I don't think you need the total, actually. The ones that are not counted by my function are the the ones  samtools skips by default (see https://www.htslib.org/doc/samtools-mpileup.html; for example, it won't count bases with a fastq quality under 13), and the ones not included in the AD value are those the variant caller filtered out. You can't just use all reads, you must only take into account the ones that pass some basic QC filtering otherwise, you just get nonsense. If you see such huge differences, I would investigate the quality (both mapping and fastq) in those regions.

Sounds good, thank you for your advice

may I ask how you produced the 1.vcf.gz file in your example?

@Nereus from one I had. I just looked at a bam file I had lying around and made the vcf point to a variant in that file and ran the command I showed you.

I see. Would it also be possible to create a custom VCF if I haven't got one available, e.g. with data from Ensembl? Thank you for your time!

@Nereus um, sure, why not? I mean, a vcf is just a text file, you can make one manually if you want.

I have noticed that your approach produces the desired result if the variant is actually present; however, when the result is not present an empty VCF file is created (I'm adding an Edit above)

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


text:
GSEA vs GSVA: Pros and cons of each method
Trying to understand when GSEA is more appropriate than GSVA and vice versa. I have seen cases when running GSEA and GSVA on the same task - compare enrichment of a geneset between two groups - gives different results. In this case one method will identify significant enrichment while the other will not. GSEA seems to detect more "hits" than GSVA.
I understand that they are different methods operating on different levels (group/condition level vs single sample level) and GSVA is needed when you want to investigate the relationship between a geneset and a non-binary grouping (other continuous variable, more than 2 categories, survival analysis etc), but I cannot find much discussion on how to interpret GSEA/GSVA results when they differ or if there are some situations for the scenario I have described (comparison of geneset enrichment in 2 groups) where one method performs better than the other. Is one method better when there are fewer samples or does one method produce more false positives than the other?
I also know there are other single sample enrichment scoring algorithms (singscore for bulk RNA-Seq, and a host of scRNA-seq methods), but I ask about GSVA as it still seems to be a popular choice in many papers.
Can someone provide justification for when to use one method over the other or point me towards any papers that explore this topic?

Do you want a literature search/pointers or some comments on why to apply one or the other?

Thoughts on when to apply on vs the other is more what I'm looking for but if this is an "unsolved" problem any papers on the the topic would be helpful. I can't seem to find much when I search around


text:
How do I summarise VCF files?
I have got a number of VCF files that all contain data on one particular variant, with every VCF being created from an individual sample. So there is one entry per file (headers omitted):
##source=Mutect2
##tumor_sample=312091_5-Ctrl1_G01_S6_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312091_5-Ctrl1_G01_S6_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=805,864|0,0;DP=1871;ECNT=1;MBQ=17,0;MFRL=273,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-3.537e+00 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:1669,0:6.794e-04:1669:252,0:233,0:1595,0:805,864,0,0

I would like to create a summary table of the parameters given within the VCFs, including of the parameters within the FORMAT field. Ideally this table would have the following as header:
Sample CHROM    POS ID  REF ALT QUAL    FILTER  GT INFO-DP FORMAT-DP AD-REF AD-ALT

and every sample/variant combination would be a line. I would like to use such a table to create statistics for my samples (on zygosity, number of ALT variants etc).
Is there a pre-made tool that is capable of doing this? I'm also open to any other suggestions. Thanks!
Edit:
Here are two VCF files (I have about 30 in total) and the desired output
First VCF:
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=FAD,Number=R,Type=Integer,Description="Count of fragments supporting each allele.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##GATKCommandLine=<ID=Mutect2,CommandLine="Mutect2 --alleles /add/2024-03-04_force_mutect/E200K.vcf.gz --output /add/2024-03-04_force_mutect/output//Ctrl5_E200K_frequency.vcf --intervals /add/2024-03-04_force_mutect/E200K.vcf.gz --input /add/2023-06-18_CJD_8_samples/312091_8-Ctrl5_H02_S3_R1_001.fastq.gz.picard.markedDup.recal.bam --reference /home/mcarta/databases/hg38.fa --f1r2-median-mq 50 --f1r2-min-bq 20 --f1r2-max-depth 200 --flow-likelihood-parallel-threads 0 --flow-likelihood-optimized-comp false --trim-to-haplotype true --exact-matching false --flow-use-t0-tag false --flow-probability-threshold 0.003 --flow-remove-non-single-base-pair-indels false --flow-remove-one-zero-probs false --flow-quantization-bins 121 --flow-fill-empty-bins-value 0.001 --flow-symmetric-indel-probs false --flow-report-insertion-or-deletion false --flow-disallow-probs-larger-than-call false --flow-lump-probs false --flow-retain-max-n-probs-base-format false --flow-probability-scaling-factor 10 --flow-order-cycle-length 4 --keep-boundary-flows false --genotype-pon-sites false --genotype-germline-sites false --af-of-alleles-not-in-resource -1.0 --mitochondria-mode false --mutect3-training-mode false --mutect3-ref-downsample 10 --mutect3-alt-downsample 20 --mutect3-non-artifact-ratio 20 --tumor-lod-to-emit 3.0 --initial-tumor-lod 2.0 --pcr-snv-qual 40 --pcr-indel-qual 40 --base-qual-correction-factor 5 --max-population-af 0.01 --downsampling-stride 1 --callable-depth 10 --max-suspicious-reads-per-alignment-start 0 --normal-lod 2.2 --ignore-itr-artifacts false --gvcf-lod-band -2.5 --gvcf-lod-band -2.0 --gvcf-lod-band -1.5 --gvcf-lod-band -1.0 --gvcf-lod-band -0.5 --gvcf-lod-band 0.0 --gvcf-lod-band 0.5 --gvcf-lod-band 1.0 --minimum-allele-fraction 0.0 --independent-mates false --flow-mode NONE --disable-adaptive-pruning false --kmer-size 10 --kmer-size 25 --dont-increase-kmer-sizes-for-cycles false --allow-non-unique-kmers-in-ref false --num-pruning-samples 1 --min-dangling-branch-length 4 --recover-all-dangling-branches false --max-num-haplotypes-in-population 128 --min-pruning 2 --adaptive-pruning-initial-error-rate 0.001 --pruning-lod-threshold 2.302585092994046 --pruning-seeding-lod-threshold 9.210340371976184 --max-unpruned-variants 100 --linked-de-bruijn-graph false --disable-artificial-haplotype-recovery false --enable-legacy-graph-cycle-detection false --debug-assembly false --debug-graph-transformations false --capture-assembly-failure-bam false --num-matching-bases-in-dangling-end-to-recover -1 --error-correction-log-odds -Infinity --error-correct-reads false --kmer-length-for-read-error-correction 25 --min-observations-for-kmer-to-be-solid 20 --likelihood-calculation-engine PairHMM --base-quality-score-threshold 18 --dragstr-het-hom-ratio 2 --dont-use-dragstr-pair-hmm-scores false --pair-hmm-gap-continuation-penalty 10 --expected-mismatch-rate-for-read-disqualification 0.02 --pair-hmm-implementation FASTEST_AVAILABLE --pcr-indel-model CONSERVATIVE --phred-scaled-global-read-mismapping-rate 45 --disable-symmetric-hmm-normalizing false --disable-cap-base-qualities-to-map-quality false --enable-dynamic-read-disqualification-for-genotyping false --dynamic-read-disqualification-threshold 1.0 --native-pair-hmm-threads 4 --native-pair-hmm-use-double-precision false --flow-hmm-engine-min-indel-adjust 6 --flow-hmm-engine-flat-insertion-penatly 45 --flow-hmm-engine-flat-deletion-penatly 45 --pileup-detection false --use-pdhmm false --use-pdhmm-overlap-optimization false --make-determined-haps-from-pd-code false --print-pileupcalling-status false --fallback-gga-if-pdhmm-fails true --pileup-detection-enable-indel-pileup-calling false --pileup-detection-active-region-phred-threshold 0.0 --num-artificial-haplotypes-to-add-per-allele 5 --artifical-haplotype-filtering-kmer-size 10 --pileup-detection-snp-alt-threshold 0.1 --pileup-detection-indel-alt-threshold 0.1 --pileup-detection-absolute-alt-depth 0.0 --pileup-detection-snp-adjacent-to-assembled-indel-range 5 --pileup-detection-snp-basequality-filter 12 --pileup-detection-bad-read-tolerance 0.0 --pileup-detection-proper-pair-read-badness true --pileup-detection-edit-distance-read-badness-threshold 0.08 --pileup-detection-chimeric-read-badness true --pileup-detection-template-mean-badness-threshold 0.0 --pileup-detection-template-std-badness-threshold 0.0 --pileup-detection-filter-assembly-alt-bad-read-tolerance 0.0 --pileup-detection-edit-distance-read-badness-for-assembly-filtering-threshold 0.12 --bam-writer-type CALLED_HAPLOTYPES --dont-use-soft-clipped-bases false --override-fragment-softclip-check false --min-base-quality-score 10 --smith-waterman FASTEST_AVAILABLE --emit-ref-confidence NONE --max-mnp-distance 1 --force-call-filtered-alleles false --reference-model-deletion-quality 30 --soft-clip-low-quality-ends false --allele-informative-reads-overlap-margin 2 --smith-waterman-dangling-end-match-value 25 --smith-waterman-dangling-end-mismatch-penalty -50 --smith-waterman-dangling-end-gap-open-penalty -110 --smith-waterman-dangling-end-gap-extend-penalty -6 --smith-waterman-haplotype-to-reference-match-value 200 --smith-waterman-haplotype-to-reference-mismatch-penalty -150 --smith-waterman-haplotype-to-reference-gap-open-penalty -260 --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 --smith-waterman-read-to-haplotype-match-value 10 --smith-waterman-read-to-haplotype-mismatch-penalty -15 --smith-waterman-read-to-haplotype-gap-open-penalty -30 --smith-waterman-read-to-haplotype-gap-extend-penalty -5 --flow-assembly-collapse-hmer-size 0 --flow-assembly-collapse-partial-mode false --flow-filter-alleles false --flow-filter-alleles-qual-threshold 30.0 --flow-filter-alleles-sor-threshold 3.0 --flow-filter-lone-alleles false --flow-filter-alleles-debug-graphs false --min-assembly-region-size 50 --max-assembly-region-size 300 --active-probability-threshold 0.002 --max-prob-propagation-distance 50 --force-active false --assembly-region-padding 100 --padding-around-indels 75 --padding-around-snps 20 --padding-around-strs 75 --max-extension-into-assembly-region-padding-legacy 25 --max-reads-per-alignment-start 50 --enable-legacy-assembly-region-trimming false --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false --max-read-length 2147483647 --min-read-length 30 --minimum-mapping-quality 20 --disable-tool-default-annotations false --enable-all-annotations false",Version="4.5.0.0",Date="March 5, 2024 at 10:20:03 AM UTC">
##INFO=<ID=AS_SB_TABLE,Number=1,Type=String,Description="Allele-specific forward/reverse read counts for strand bias tests. Includes the reference and alleles separated by |.">
##INFO=<ID=AS_UNIQ_ALT_READ_COUNT,Number=A,Type=Integer,Description="Number of reads with unique start and mate end positions for each alt at a variant site">
##INFO=<ID=CONTQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to contamination">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=ECNT,Number=1,Type=Integer,Description="Number of events in this haplotype">
##INFO=<ID=GERMQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not germline variants">
##INFO=<ID=MBQ,Number=R,Type=Integer,Description="median base quality by allele">
##INFO=<ID=MFRL,Number=R,Type=Integer,Description="median fragment length by allele">
##INFO=<ID=MMQ,Number=R,Type=Integer,Description="median mapping quality by allele">
##INFO=<ID=MPOS,Number=A,Type=Integer,Description="median distance from end of read">
##INFO=<ID=NALOD,Number=A,Type=Float,Description="Negative log 10 odds of artifact in normal with same allele fraction as tumor">
##INFO=<ID=NCount,Number=1,Type=Integer,Description="Count of N bases in the pileup">
##INFO=<ID=NLOD,Number=A,Type=Float,Description="Normal log 10 likelihood ratio of diploid het or hom alt genotypes">
##INFO=<ID=OCM,Number=1,Type=Integer,Description="Number of alt reads whose original alignment doesn't match the current contig.">
##INFO=<ID=PON,Number=0,Type=Flag,Description="site found in panel of normals">
##INFO=<ID=POPAF,Number=A,Type=Float,Description="negative log 10 population allele frequencies of alt alleles">
##INFO=<ID=ROQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to read orientation artifact">
##INFO=<ID=RPA,Number=R,Type=Integer,Description="Number of times tandem repeat unit is repeated, for each allele (including reference)">
##INFO=<ID=RU,Number=1,Type=String,Description="Tandem repeat unit (bases)">
##INFO=<ID=SEQQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not sequencing errors">
##INFO=<ID=STR,Number=0,Type=Flag,Description="Variant is a short tandem repeat">
##INFO=<ID=STRANDQ,Number=1,Type=Integer,Description="Phred-scaled quality of strand bias artifact">
##INFO=<ID=STRQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles in STRs are not polymerase slippage errors">
##INFO=<ID=TLOD,Number=A,Type=Float,Description="Log 10 likelihood ratio score of variant existing versus not existing">
##MutectVersion=2.2
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr4,length=190214555>
##filtering_status=Warning: unfiltered Mutect 2 calls.  Please run FilterMutectCalls to remove false positives.
##source=Mutect2
##tumor_sample=312091_8-Ctrl5_H02_S3_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312091_8-Ctrl5_H02_S3_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=1376,1402|0,0;DP=2939;ECNT=1;MBQ=18,0;MFRL=257,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-3.455e+00   GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:2778,0:4.100e-04:2778:472,0:398,0:2493,0:1376,1402,0,0

Second VCF:
##fileformat=VCFv4.2
##FORMAT=<ID=AD,Number=R,Type=Integer,Description="Allelic depths for the ref and alt alleles in the order listed">
##FORMAT=<ID=AF,Number=A,Type=Float,Description="Allele fractions of alternate alleles in the tumor">
##FORMAT=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth (reads with MQ=255 or with bad mates are filtered)">
##FORMAT=<ID=F1R2,Number=R,Type=Integer,Description="Count of reads in F1R2 pair orientation supporting each allele">
##FORMAT=<ID=F2R1,Number=R,Type=Integer,Description="Count of reads in F2R1 pair orientation supporting each allele">
##FORMAT=<ID=FAD,Number=R,Type=Integer,Description="Count of fragments supporting each allele.">
##FORMAT=<ID=GQ,Number=1,Type=Integer,Description="Genotype Quality">
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
##FORMAT=<ID=PGT,Number=1,Type=String,Description="Physical phasing haplotype information, describing how the alternate alleles are phased in relation to one another; will always be heterozygous and is not intended to describe called alleles">
##FORMAT=<ID=PID,Number=1,Type=String,Description="Physical phasing ID information, where each unique ID within a given sample (but not across samples) connects records within a phasing group">
##FORMAT=<ID=PL,Number=G,Type=Integer,Description="Normalized, Phred-scaled likelihoods for genotypes as defined in the VCF specification">
##FORMAT=<ID=PS,Number=1,Type=Integer,Description="Phasing set (typically the position of the first variant in the set)">
##FORMAT=<ID=SB,Number=4,Type=Integer,Description="Per-sample component statistics which comprise the Fisher's Exact Test to detect strand bias.">
##GATKCommandLine=<ID=Mutect2,CommandLine="Mutect2 --alleles /add/2024-03-04_force_mutect/E200K.vcf.gz --output /add/2024-03-04_force_mutect/output//Ctrl4_E200K_frequency.vcf --intervals /add/2024-03-04_force_mutect/E200K.vcf.gz --input /add/2023-06-18_CJD_16_samples/312092_15-Ctrl4_H01_S11_R1_001.fastq.gz.picard.markedDup.recal.bam --reference /home/mcarta/databases/hg38.fa --f1r2-median-mq 50 --f1r2-min-bq 20 --f1r2-max-depth 200 --flow-likelihood-parallel-threads 0 --flow-likelihood-optimized-comp false --trim-to-haplotype true --exact-matching false --flow-use-t0-tag false --flow-probability-threshold 0.003 --flow-remove-non-single-base-pair-indels false --flow-remove-one-zero-probs false --flow-quantization-bins 121 --flow-fill-empty-bins-value 0.001 --flow-symmetric-indel-probs false --flow-report-insertion-or-deletion false --flow-disallow-probs-larger-than-call false --flow-lump-probs false --flow-retain-max-n-probs-base-format false --flow-probability-scaling-factor 10 --flow-order-cycle-length 4 --keep-boundary-flows false --genotype-pon-sites false --genotype-germline-sites false --af-of-alleles-not-in-resource -1.0 --mitochondria-mode false --mutect3-training-mode false --mutect3-ref-downsample 10 --mutect3-alt-downsample 20 --mutect3-non-artifact-ratio 20 --tumor-lod-to-emit 3.0 --initial-tumor-lod 2.0 --pcr-snv-qual 40 --pcr-indel-qual 40 --base-qual-correction-factor 5 --max-population-af 0.01 --downsampling-stride 1 --callable-depth 10 --max-suspicious-reads-per-alignment-start 0 --normal-lod 2.2 --ignore-itr-artifacts false --gvcf-lod-band -2.5 --gvcf-lod-band -2.0 --gvcf-lod-band -1.5 --gvcf-lod-band -1.0 --gvcf-lod-band -0.5 --gvcf-lod-band 0.0 --gvcf-lod-band 0.5 --gvcf-lod-band 1.0 --minimum-allele-fraction 0.0 --independent-mates false --flow-mode NONE --disable-adaptive-pruning false --kmer-size 10 --kmer-size 25 --dont-increase-kmer-sizes-for-cycles false --allow-non-unique-kmers-in-ref false --num-pruning-samples 1 --min-dangling-branch-length 4 --recover-all-dangling-branches false --max-num-haplotypes-in-population 128 --min-pruning 2 --adaptive-pruning-initial-error-rate 0.001 --pruning-lod-threshold 2.302585092994046 --pruning-seeding-lod-threshold 9.210340371976184 --max-unpruned-variants 100 --linked-de-bruijn-graph false --disable-artificial-haplotype-recovery false --enable-legacy-graph-cycle-detection false --debug-assembly false --debug-graph-transformations false --capture-assembly-failure-bam false --num-matching-bases-in-dangling-end-to-recover -1 --error-correction-log-odds -Infinity --error-correct-reads false --kmer-length-for-read-error-correction 25 --min-observations-for-kmer-to-be-solid 20 --likelihood-calculation-engine PairHMM --base-quality-score-threshold 18 --dragstr-het-hom-ratio 2 --dont-use-dragstr-pair-hmm-scores false --pair-hmm-gap-continuation-penalty 10 --expected-mismatch-rate-for-read-disqualification 0.02 --pair-hmm-implementation FASTEST_AVAILABLE --pcr-indel-model CONSERVATIVE --phred-scaled-global-read-mismapping-rate 45 --disable-symmetric-hmm-normalizing false --disable-cap-base-qualities-to-map-quality false --enable-dynamic-read-disqualification-for-genotyping false --dynamic-read-disqualification-threshold 1.0 --native-pair-hmm-threads 4 --native-pair-hmm-use-double-precision false --flow-hmm-engine-min-indel-adjust 6 --flow-hmm-engine-flat-insertion-penatly 45 --flow-hmm-engine-flat-deletion-penatly 45 --pileup-detection false --use-pdhmm false --use-pdhmm-overlap-optimization false --make-determined-haps-from-pd-code false --print-pileupcalling-status false --fallback-gga-if-pdhmm-fails true --pileup-detection-enable-indel-pileup-calling false --pileup-detection-active-region-phred-threshold 0.0 --num-artificial-haplotypes-to-add-per-allele 5 --artifical-haplotype-filtering-kmer-size 10 --pileup-detection-snp-alt-threshold 0.1 --pileup-detection-indel-alt-threshold 0.1 --pileup-detection-absolute-alt-depth 0.0 --pileup-detection-snp-adjacent-to-assembled-indel-range 5 --pileup-detection-snp-basequality-filter 12 --pileup-detection-bad-read-tolerance 0.0 --pileup-detection-proper-pair-read-badness true --pileup-detection-edit-distance-read-badness-threshold 0.08 --pileup-detection-chimeric-read-badness true --pileup-detection-template-mean-badness-threshold 0.0 --pileup-detection-template-std-badness-threshold 0.0 --pileup-detection-filter-assembly-alt-bad-read-tolerance 0.0 --pileup-detection-edit-distance-read-badness-for-assembly-filtering-threshold 0.12 --bam-writer-type CALLED_HAPLOTYPES --dont-use-soft-clipped-bases false --override-fragment-softclip-check false --min-base-quality-score 10 --smith-waterman FASTEST_AVAILABLE --emit-ref-confidence NONE --max-mnp-distance 1 --force-call-filtered-alleles false --reference-model-deletion-quality 30 --soft-clip-low-quality-ends false --allele-informative-reads-overlap-margin 2 --smith-waterman-dangling-end-match-value 25 --smith-waterman-dangling-end-mismatch-penalty -50 --smith-waterman-dangling-end-gap-open-penalty -110 --smith-waterman-dangling-end-gap-extend-penalty -6 --smith-waterman-haplotype-to-reference-match-value 200 --smith-waterman-haplotype-to-reference-mismatch-penalty -150 --smith-waterman-haplotype-to-reference-gap-open-penalty -260 --smith-waterman-haplotype-to-reference-gap-extend-penalty -11 --smith-waterman-read-to-haplotype-match-value 10 --smith-waterman-read-to-haplotype-mismatch-penalty -15 --smith-waterman-read-to-haplotype-gap-open-penalty -30 --smith-waterman-read-to-haplotype-gap-extend-penalty -5 --flow-assembly-collapse-hmer-size 0 --flow-assembly-collapse-partial-mode false --flow-filter-alleles false --flow-filter-alleles-qual-threshold 30.0 --flow-filter-alleles-sor-threshold 3.0 --flow-filter-lone-alleles false --flow-filter-alleles-debug-graphs false --min-assembly-region-size 50 --max-assembly-region-size 300 --active-probability-threshold 0.002 --max-prob-propagation-distance 50 --force-active false --assembly-region-padding 100 --padding-around-indels 75 --padding-around-snps 20 --padding-around-strs 75 --max-extension-into-assembly-region-padding-legacy 25 --max-reads-per-alignment-start 50 --enable-legacy-assembly-region-trimming false --interval-set-rule UNION --interval-padding 0 --interval-exclusion-padding 0 --interval-merging-rule ALL --read-validation-stringency SILENT --seconds-between-progress-updates 10.0 --disable-sequence-dictionary-validation false --create-output-bam-index true --create-output-bam-md5 false --create-output-variant-index true --create-output-variant-md5 false --max-variants-per-shard 0 --lenient false --add-output-sam-program-record true --add-output-vcf-command-line true --cloud-prefetch-buffer 40 --cloud-index-prefetch-buffer -1 --disable-bam-index-caching false --sites-only-vcf-output false --help false --version false --showHidden false --verbosity INFO --QUIET false --use-jdk-deflater false --use-jdk-inflater false --gcs-max-retries 20 --gcs-project-for-requester-pays  --disable-tool-default-read-filters false --max-read-length 2147483647 --min-read-length 30 --minimum-mapping-quality 20 --disable-tool-default-annotations false --enable-all-annotations false",Version="4.5.0.0",Date="March 5, 2024 at 10:21:02 AM UTC">
##INFO=<ID=AS_SB_TABLE,Number=1,Type=String,Description="Allele-specific forward/reverse read counts for strand bias tests. Includes the reference and alleles separated by |.">
##INFO=<ID=AS_UNIQ_ALT_READ_COUNT,Number=A,Type=Integer,Description="Number of reads with unique start and mate end positions for each alt at a variant site">
##INFO=<ID=CONTQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to contamination">
##INFO=<ID=DP,Number=1,Type=Integer,Description="Approximate read depth; some reads may have been filtered">
##INFO=<ID=ECNT,Number=1,Type=Integer,Description="Number of events in this haplotype">
##INFO=<ID=GERMQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not germline variants">
##INFO=<ID=MBQ,Number=R,Type=Integer,Description="median base quality by allele">
##INFO=<ID=MFRL,Number=R,Type=Integer,Description="median fragment length by allele">
##INFO=<ID=MMQ,Number=R,Type=Integer,Description="median mapping quality by allele">
##INFO=<ID=MPOS,Number=A,Type=Integer,Description="median distance from end of read">
##INFO=<ID=NALOD,Number=A,Type=Float,Description="Negative log 10 odds of artifact in normal with same allele fraction as tumor">
##INFO=<ID=NCount,Number=1,Type=Integer,Description="Count of N bases in the pileup">
##INFO=<ID=NLOD,Number=A,Type=Float,Description="Normal log 10 likelihood ratio of diploid het or hom alt genotypes">
##INFO=<ID=OCM,Number=1,Type=Integer,Description="Number of alt reads whose original alignment doesn't match the current contig.">
##INFO=<ID=PON,Number=0,Type=Flag,Description="site found in panel of normals">
##INFO=<ID=POPAF,Number=A,Type=Float,Description="negative log 10 population allele frequencies of alt alleles">
##INFO=<ID=ROQ,Number=1,Type=Float,Description="Phred-scaled qualities that alt allele are not due to read orientation artifact">
##INFO=<ID=RPA,Number=R,Type=Integer,Description="Number of times tandem repeat unit is repeated, for each allele (including reference)">
##INFO=<ID=RU,Number=1,Type=String,Description="Tandem repeat unit (bases)">
##INFO=<ID=SEQQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles are not sequencing errors">
##INFO=<ID=STR,Number=0,Type=Flag,Description="Variant is a short tandem repeat">
##INFO=<ID=STRANDQ,Number=1,Type=Integer,Description="Phred-scaled quality of strand bias artifact">
##INFO=<ID=STRQ,Number=1,Type=Integer,Description="Phred-scaled quality that alt alleles in STRs are not polymerase slippage errors">
##INFO=<ID=TLOD,Number=A,Type=Float,Description="Log 10 likelihood ratio score of variant existing versus not existing">
##MutectVersion=2.2
##contig=<ID=chr2,length=242193529>
##contig=<ID=chr20,length=64444167>
##contig=<ID=chr4,length=190214555>
##filtering_status=Warning: unfiltered Mutect 2 calls.  Please run FilterMutectCalls to remove false positives.
##source=Mutect2
##tumor_sample=312092_15-Ctrl4_H01_S11_R1_001.fastq.gz
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  312092_15-Ctrl4_H01_S11_R1_001.fastq.gz
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=22,12|0,0;DP=48;ECNT=1;MBQ=11,0;MFRL=243,0;MMQ=60,60;MPOS=50;POPAF=7.30;TLOD=-1.584e+00 GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:34,0:0.026:34:0,0:0,0:33,0:22,12,0,0

Desired output (values from the first VCF above, with tabs as delimiter)
Sample CHROM    POS ID  REF ALT QUAL    FILTER  GT INFO-DP FORMAT-DP AD-REF AD-ALT
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz chr20 4699818 . G A . . 0/1 2939 2778 2778 0

I suppose merging them would be an option. Is there a method to extract the information I'd like from a merged VCF?

The first step is to merge your VCFs. For this, you need to compress them with bgzip and then index with tabix:
for file in *vcf; do
  bgzip "$file"
  tabix "$file".gz
done

Now, feed the files to bcftools merge. In this example, I have saved your two files as sample1.vcf.gz and sample2.vcf.gz, and they're the only vcf files in the directory, so I can catch them all with *vcf.gz:
bcftools merge *.vcf.gz > merged.vcf

This results in a file containing a single variant line, with an entry for each of the two samples:
chr20   4699818 .   G   A   .   .   AS_SB_TABLE=1376,1402|0,0;ECNT=1;DP=2987;MBQ=11,0;MFRL=243,0;MMQ=60,60;MPOS=50;POPAF=7.3;TLOD=-1.584    GT:AD:AF:DP:F1R2:F2R1:FAD:SB    0/1:2778,0:0.00041:2778:472,0:398,0:2493,0:1376,1402,0,0    0/1:34,0:0.026:34:0,0:0,0:33,0:22,12,0,0

From this, bcftools query can give you almost exactly what you need. The main problem is that when merging, the INFO fields are not combined. Only one is kept. In my test, it was the one from the last file processed, but this could be random, I don't know. This is the closest I could get to your desired format:
$ { printf 'Sample\tCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tGT\tFORMAT-DP\tAD\n';
  bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf ; } | column -s$'\t' -t
Sample                                                                          CHROM  POS      ID  REF  ALT  QUAL  FILTER  GT          FORMAT-DP  AD
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz 312092_15-Ctrl4_H01_S11_R1_001.fastq.gz   chr20  4699817  .   G    A    .     .        0/1  0/1   2778 34    2778,0 34,0

The main command is this:
bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf

The rest just add a header and pretty formatting for clarity. Please see https://samtools.github.io/bcftools/bcftools.html#query for details on the options used.
This is close, but not quite. You have the sample names listed on the same line, we're missing the DP value from the FORMAT field and we have the ADs shown in pairs, not in different columns as requested. A little post-processing can fix everything except the missing DP/FORMAT:
$ {
  printf 'Sample\tCHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tGT\tFORMAT-DP\tAD-REF\tAD-ALT\n';
  bcftools query --format "[%SAMPLE ]\t%CHROM\t%POS0\t%ID\t%REF\t%ALT\t%QUAL\t%FILTER\t[%GT ]\t[%DP ]\t[%AD ] " merged.vcf |
   awk -F'\t' -v OFS="\t" '
  {
    split($1, samples, " ");
    split($9, genotypes, " ");
    split($10, dps, " ");
    split($11, ads, " ");
    for(s in samples){
       split(ads[s], sampleAd, ",");
       print samples[s], $2, $3, $4, $5, $6, $7, $8, genotypes[s], dps[s], sampleAd[1], sampleAd[2]
     }
  }';
} | column -s$'\t' -t
Sample                                   CHROM  POS      ID       REF  ALT  QUAL  FILTER  GT   FORMAT-DP  AD-REF  AD-ALT
312091_8-Ctrl5_H02_S3_R1_001.fastq.gz    chr20  4699817  .  G    A    .     .       0/1  2778       2778    0
312092_15-Ctrl4_H01_S11_R1_001.fastq.gz  chr20  4699817 .  G    A    .     .       0/1  34         34      0

Important: the command above assumes there are no spaces in your sample names nor in any of the values returned. It also assumes only one ALT allele per variant. Please update your question if these assumptions are not valid.
This is as close as I can get it using "standard" tools and relatively little post-processing. It is everything you asked for except the INFO DP. If you absolutely must have that too (which seems a bit odd, honestly, you usually want the one actually reported for the variant), I fear that will require just parsing the whole thing manually.

this is exactly what I needed, thank you so much!

The simplest alternative is to have a directory with all the vcf files you want to merge, and use bcftools merge. Example:
vcf_dir=path_to_your_vcf_directory
bcftools merge $vcf_dir/*.vcf -o bcftools_merged.vcf

This does not answer the question. It's literally the second step from terdon's answer, so I don't see why this deserves to be its own answer.

Just to point out @RamRS does have a valid complaint.


text:
How to extract mutations in a particular gene from the VCF file of whole genome
Also posted on biostars
I would like to extract all the mutations in a single gene of interest from the VCF file. For example I would like to extract all the mutations in the NFKB1 gene from the whole genome VCF file.

What have you tried? bcftools view has an option to extract by coordinates or by annotation if your VCF had the gene annotation. If you used custom annotation tools, they may have their own filtering utilities.

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.

Find the location of the gene of interest and then filter the vcf for that location using vcftools, bcftools, plink, or some other software. There is a lot of good documentation available online for these tools.

As mentioned, there are many options to go (eg. vcftools and bcftools). The Broad Institute provides a valuable toolset (GATK) for this very purpose and extended procedures for various variant selection tasks.
SelectVariants function can be good for this purpose, though the location of a gene of interest should be addressed:
https://gatk.broadinstitute.org/hc/en-us/articles/13832694334235-SelectVariants


text:
Why is mtDNA not labelled as such for many organisms?
Also posted on biostars
While I was searching the genomes of some eukariotic organisms (fungi,plants,protists) even though I knew that they had mitochondrial DNA,when I searched their genome structure I couldn't find any "chromosome" labeled as "mtDNA". Here I give some links of genomes that show this:
https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF_016861865.1
https://www.ncbi.nlm.nih.gov/genome/gdv/browser/genome/?id=GCF_009556855.1
Why does this happens? How can I find if the mtDNA sequence is inluded in these "chromosomes"?

Note how in https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_016861865.1/, the number of "organelles" is listed as 0 for RefSeq but 1 in GenBank. I would guess you can find a mitochondrial sequence in GenBank.

The most obvious answer is wet-lab rather than bioinformatics: mtDNA can get removed in the DNA extraction process. As a control group, in bacterial drug-resistance studies the plasmids are frequently lost. mtDNA at ~15KB is smaller than a regular bacterial plasmid.
On the bioinformatics front there would be an issue about reference genome assemblies as well, i.e. if mtDNA genome isn't present then it will not appear as a separate genome.


text:
All against all sequence alignment
I have a large set of protein sequences and would like to create a "similarity matrix" out of them (for phylo trees later). So I have to run n(n-1) sequence alignments. Wondering what library I should use for that and also what sophisticated lib would provide this all against all without me writing a for loop and parsing results myself...
Currently using Bio.pairwise2.align.globalxx from BioPython as a development placeholder but I do not think the results look really useful? It returns a score > 1 and no similarity in [0,1[ or evalue.
I have been out of the Bioinfo game for a while, so any suggestions for libraries (python preferred) much appreciated! :)

I wouldn't use Python here. I've looked at using Bio.pairwise2.align. and its has major shortcomings. In this context the problem is establishing homologous alignment positions ... if you want to make a tree.
Whats the problems? This is you start introducing the gap character '-' as an additional state in constructing your tree and the tree will bias towards not the maximum number of substitutions but the maximum gap (indel) difference between the sequences.
The second problem is changes in homologous sites between different pairwise alignments. This gets complicated fast.
Is that bad? The problem is we have little understanding of the mutation model behind generating indels and there isn't a model to account for varying descriptions of homologous alignment positions. However, we have a very clear understanding of amino acid substitutions in a single global alignment, i.e. a singular description of homologous amino acid positions. A lot of work has gone into modelling amino acid mutation, with a tree a practical output of that work.

Thus, I would simply align using muscle5 under the -super5 option. This will produce a global alignment very quickly and for vast amounts of data. This alignment will produce the best model of homologous amino acid sites.
If you just want a distance matrix then import the entire alignment into MEGAXI and request the JTT model of amino acid substitution with the gamma distribution. This will output a distance matrix nicely and will then allow you to produce a tree.
BTW estimating the gamma distribution without maximum likelihood can be difficult in distance methods (MEGA might have resolved this). This tree will be much better than the one you are proposing (you can estimate this parameter via e.g. IQTREE or ModelTest).
Why? Again alot of work has been done on modelling amino acid substations and JTT is one of the better models (assuming this is nuclear sequence not mtDNA). The LG matrix is much better BTW, but I don't think thats in MEGA.


text:
Tools and software to visualize the moving cell compartments in 3D
Also asked on biostars
I am looking for tools and software that could help me visualize moving cell compartments in 3D where I can go from the cell membrane to the ER, then out with the protein into the matrix. In short a 3D visualization of the cell environment.
Since I am a newbie, I would appreciate any tool/software suggestion for me to try.

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.


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.


text:
Exploring Amino Acid Patterns in Proteins Through N-gram Analysis
In our recent research, we delved into the amino acid composition of protein sequences by applying n-gram analysis techniques. By examining the distribution of n-grams ranging from 1-gram to 11-gram, we aimed to discern underlying patterns that could be indicative of structural or functional significance.
We started with an exploration of 1-grams to understand the basic composition of our protein dataset, followed by investigating the prevalence and distribution of larger n-grams. Our analysis included evaluating the fit of these distributions against well-known statistical laws—namely, Benford's, Pareto, and Zipf's.
An intriguing pattern emerged as n-gram lengths increased. The frequency distributions began to approximate Benford's Law more closely, particularly at tetragram and pentagram levels. Beyond pentagrams, the distributions deviated, suggesting a higher level of sequence diversity and complexity.
Furthermore, we assessed how these patterns aligned with the Pareto principle. While we found that the data did not strictly adhere to the "80/20 rule", there was an interesting variation in the concentration of occurrences across different n-gram lengths.
Our approach also involved identifying sequences with anomalies such as the presence of consecutive 'X' characters, which denotes unknown or unspecified amino acids, and sequences that were unusually short or long compared to the general protein population.
As we continue to dissect these patterns, several questions arise, inviting further scrutiny and discussion within the scientific community:

What biological insights can be inferred from the observed fit of amino acid patterns to Benford's Law, particularly for tetragrams and pentagrams?
How might the deviations in longer n-grams inform our understanding of protein complexity and functionality?
In what ways can the prevalence of specific n-grams guide the development of more accurate predictive models for protein function?
What are the implications of the identified sequence anomalies for data quality and sequence annotation in protein databases?

We are eager to engage with the bioinformatics community to explore these questions and welcome any insights or collaborative ideas that can drive this research forward.
Our public notebook Testing Pareto, Benford and Zipf

This question currently includes multiple questions in one. It should focus on one specific problem only. Note that Bioinformatics Stack Exchange posts are not a discussion forum; try chat or reddit for that.


text:
Creating custom VCF file
I'm stuck on what seems to be a banal problem, but I can't work out where the issue is. I need to create some customised VCF files that specify particular variants of interest.
This is a custom VCF I created for this variant:
##fileformat=VCFv4.3
##FILTER=<ID=PASS,Description="All filters passed">
##fileDate=20240308
##source=Ensembl
##reference=GRCh38
#CHROM  POS ID  REF ALT QUAL  FILTER    INFO
20  4699752 rs74315403  G   A   .   .  .

The formatting seems correct to me; the values are separated by tabs, all relevant column headers are present, there is a value for every header. I try to use it in the Mutect2 --alleles and -L options, as follows:
gatk Mutect2 \
  -R ref.fa \
  --alleles my.vcf.gz  \
  -L my.vcf.gz\
  -I input.bam \
  -O output.vcf

This gives me the following error:
Caused by: htsjdk.tribble.TribbleException$MalformedFeatureFile: Unable to parse header with error: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM     POS     ID      REF     ALT     QUAL  FILTER    INFO, for input source: path/to/my.vcf
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:265)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:104)
        at htsjdk.tribble.TribbleIndexedFeatureReader.<init>(TribbleIndexedFeatureReader.java:129)
        at htsjdk.tribble.AbstractFeatureReader.getFeatureReader(AbstractFeatureReader.java:121)
        at org.broadinstitute.hellbender.engine.FeatureDataSource.getTribbleFeatureReader(FeatureDataSource.java:433)
        ... 15 more
Caused by: htsjdk.tribble.TribbleException$InvalidHeader: Your input file has a malformed header: there are not enough columns present in the header line: #CHROM       POS     ID      REF    ALT      QUAL  FILTER    INFO
        at htsjdk.variant.vcf.AbstractVCFCodec.parseHeaderFromLines(AbstractVCFCodec.java:149)
        at htsjdk.variant.vcf.VCFCodec.readActualHeader(VCFCodec.java:117)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:79)
        at htsjdk.tribble.AsciiFeatureCodec.readHeader(AsciiFeatureCodec.java:37)
        at htsjdk.tribble.TribbleIndexedFeatureReader.readHeader(TribbleIndexedFeatureReader.java:263)

Supposedly there's something wrong with the header. Can you spot the issue? Thanks!

That sounds like your header isn't tab separated. You can check for the presence of tabs with od -c. So, to check the header, do
zgrep '^#C' my.vcf.gz | od -c

Then check the output, it should show \t between each header field, like this:
$ zgrep '^#C'  variants.vcf.gz | od -c
0000000   #   C   H   R   O   M  \t   P   O   S  \t   I   D  \t   R   E
0000020   F  \t   A   L   T  \t   Q   U   A   L  \t   F   I   L   T   E
0000040   R  \t   I   N   F   O  \t   F   O   R   M   A   T  \t   E   D
0000060   2   0   0   1   3   9   6  \n
0000070

Nice on the od -c. I generally use cat -te, but this looks easier to eyeball


text:
Where would I get the activity profile dataset for Lantus Solostar insulin, and for other diabetes treatments such as Metformin?
I am working on a personal project for a better, more informed management of my own Diabetes Mellitus, Type 2 (aka "T2DM", or "T2D"). With some difficulty I am able to access and parse my Interstitial Glucose Level ("IGL") from my Continuous Glucose Monitor ("CGM"), and import the minute-by-minute data into a spreadsheet. I would like to be able to compare my IGL results to the activity profiles of both my prescribed insulin and my prescribed oral diabetes medications. I anticipate that such a comparison of the IGL and activity profiles of the treatments should give me an insight into the best times for administering my treatments. Later I will be adding further layers of comparison to this, such as my blood pressure, pulse, blood oxygen SpO₂ level, carbohydrate intake, exercise, etc to get a statistical view of which of these factors influence the IGL numbers the most, over time.
What I need first to be able to do all of this is a copy of the dataset for the activity profile of the insulin that I have been prescribed, the Lantus SoloStar 100ɪᴜ/mL pen. I have seen the graph of this data available online quite often, but I need to get the actual figures that the graphs are based on, in order to use them in my own modelling. If anyone knows where I might obtain this information, as a non-clinician myself, please let me know. Similarly, for a source of such activity profile data for common diabetes drugs such as Metformin and Gliclazide.
Also, if anyone knows of a better way to parse this information, including the raw data from the Abbott Freestyle Libre 2 with LibreView, please do let me know. Thanks.

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.

You mention that you have seen the graph of the data you need online but you need the actual data.  Can you provide a link to the example graph? Ideally you would find a scientific article that shows that graph, and then look in the supplemental material for a table used to create the graph.


text:
Compare FASTQ reconstructions to mtDNA
I have some doubts in how to proceed to compare ancient mtDNA with modern mtDNA. A little more context:
The problem:
I have some FASTQs associated with a given FASTA. They contain mtDNA. My goal is to extract the mtDNA from those FASTQs and compare with a reference mtDNA to calculate the average and SD of that difference.
I currently download the samples FASTA and FASTQs from AMTDB (https://amtdb.org/samples).
My final goal is to compare ancient mtDNA with modern mtDNA to analyze the deterioration based on the degradation patterns.
What I'm lacking in knowledge:

Does it make sense if I align the FASTQ with the corresponded downloaded FASTA? I mean, if the FASTA was already generated with those reads, wouldn't that be redundant?

Could I use an arbitrary mtDNA FASTA to align and compare? That would assume mtDNA in humans is very similar (I have doubts if that would be a safe assumption, since my goal is to compare those ancient samples with a modern and not damaged sample)

In any case, is there a recommended program to compare the results of the processing of the FASTQs with the reference FASTA?

Thank you for your help.

Can you please [edit] your question and explain the relationship between these fastqs and fasta files you mention? You should probably also tell us what species you are working with, should we assume human?

Please edit the question to limit it to a specific problem with enough detail to identify an adequate answer.


text:
How to interpret high support values for terminal nodes and low for more ancient nodes?
I made a phylogenetic tree for about 300 species using a region of the COX1 gene. The tree looks visually meaningful: species of the same genus are all clustered together.
Then I performed a 350-replicate bootstrapping and I analyzed the support values. Surprisingly I found very low values for nodes that are more close to the root and high support values for nodes that are more close to the leaves.
Almost all nodes that split into 2 leaves have Bootstrap >95% but the ones close to the root have Bootstrap<20%.
Does that mean I chose a highly variable DNA region?
How would you re-build the tree to make it more consistent on more ancient nodes? I was thinking of using a concatenation of sequences with both high and low variability regions, but i am not sure it would work.


text:
Downloading multiple SRA files as FASTQ files using SRA Toolkit
I was able to download all of the SRA files for a specific bioproject using the command:
prefetch -O output_directory --option-file SRR_Acc_List.txt

I then tried the following command to download those SRA files into FASTQ files:
find -name \*.sra -exec bsub fastq-dump --split-3 --gzip {} \;

It gave the error:
find: illegal option -- n
usage: find [-H | -L | -P] [-EXdsx] [-f path] path ... [expression]
       find [-H | -L | -P] [-EXdsx] -f path [path ...] [expression]

I then tried to add the path:
find /Users/lucykelly/Desktop/Assel_paper/output_directory -name '*.sra' -exec bsub fastq-dump --split-3 --gzip {} \;

It gave the error:
find: bsub: No such file or directory

I then tried
find . -name '*.sra' -exec fasterq-dump --gzip {} \;

and it gave the error:
find: fasterq-dump: No such file or directory

I agree with @ATPoint. At first I thought it was bsub within the Slurm cluster. If you use which fasterq-dump it should throw a blank. The bsub error is weird because you should be on the cluster right? Anyway once diagnosed its easy to solve.

That is a classical problem where the folder containing the binary of the tool is not in PATH, see for example https://askubuntu.com/questions/109381/how-to-add-path-of-a-program-to-path-environment-variable or many other threads towards PATH in the SE network.

Thank you so much @ATpoint


text:
Finding insertion in a plasmid
This question was also asked on Biostars
I received a set of Nanopore long reads of a ~3kbp long known plasmid with an unknown insert between 1kbp and 7kbp. Could you suggest a pipeline to assemble the reads and find the insertion?

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.

When I've needed to do this, I use canu to assemble the nanopore reads (choosing the highest-quality reads that are at least half the length of the target), then map the assembled consensus back to the reference sequence.
The Canu assembly part of that is described here:
https://www.protocols.io/view/plasmid-sequence-analysis-from-long-reads-36wgq4n5yvk5?step=12
There are numerous different ways to map the consensus assembled sequence to the reference sequence. BLASTn in "align 2 or more sequences" mode could be used, for example. Another option, because the sequences are relatively small and you're just aligning 2, is needle from the EMBOSS suite.


text:
Protein conservation analysis with trident scores in R script
I am currently working on analyzing protein conservation across amino acid positions using Trident scores as output from the MstatX tool. I have 40,000 Trident score files, each containing amino acid positions and their respective scores.
To streamline my analysis, I initially wrote a bash script to loop through each file and then implemented an R script for the calculations. The R script calculates moving averages with different spans to identify regions with high deviation in local protein conservation, referred to as "hotspots." These hotspots are then analyzed to determine the strongest deviation, its position, and related metrics.
Here is a snippet of my R script:
# assign input files to a vector
args = commandArgs(trailingOnly = T)

# put sample name in variable
sample<-sprintf("%s", args[1])

# load sample and make data frame
raw.entropy<-read.table(sprintf("%s", args[1]))
col.names<-c("AA.pos", "Entropy.Score")
colnames(raw.entropy)<-col.names
AA.pos<-raw.entropy$AA.pos
Entropy.Score<-raw.entropy$Entropy.Score

# make sample name nicer
sub("_entropyscores", "", sample)->sample

# load ggplot
library(ggplot2)

# calculate short and long moving averages
mydata<-ggplot(data=raw.entropy, aes(x=AA.pos)) +
  geom_point(aes(y=Entropy.Score), size=2, color="black") +
  geom_smooth(mapping=aes(y=Entropy.Score), method = "loess", color = "red", span = 0.1) +
  geom_smooth(mapping=aes(y=Entropy.Score), method = "loess", color = "blue", span = 0.4) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5, face = "bold")) +
  labs(title = sprintf("%s", args[1]), x = "Amino Acid Positions", y = "Trident Score")

# infer ratios between short and long span moving averages
# ratios represent extent to which local conservation is adhered to
# ratio around 1 means perfect conservation, ratio around 0 means poor conservation
rolling_mean_10<-layer_data(mydata, 2)$y
rolling_mean_50<-layer_data(mydata, 3)$y
ratios<-rolling_mean_10/rolling_mean_50

# find index value of 10th percentile
lowerbound<-NULL
for (i in (0.1*length(rolling_mean_10)):length(rolling_mean_10)){
if (rolling_mean_10[i] >= 0.9*median(raw.entropy$Entropy.Score)){
lowerbound<-i
break
}
}

# find index value of 90th percentile
upperbound<-NULL
for (i in (0.9*length(rolling_mean_10)):1){
if (rolling_mean_10[i] >= 0.9*median(raw.entropy$Entropy.Score)){
upperbound<-i
break
}
}

# find location in sequence with highest deviation in local protein conservation
which(ratios==min(ratios[lowerbound:upperbound]))->hotspot
ratios[hotspot]->cutoff
layer_data(mydata, 3)$x[hotspot]->x
percentile<-x/length(raw.entropy$AA.pos)*100

# write relevant items to standard output
# order is cluster name, strongest deviation, rounded amino acid position, percentile where deviation occurred
cat(sprintf("%s %s %s %s %s\n", sample, cutoff, round(x), round(percentile, 2)), file="./output_deviation.txt", append=TRUE)

After writing the output file, I aim to create a density plot with the x-axis representing the strongest deviation of each input file. However, the current method does not capture the expected bimodal distribution. I am seeking advice on alternative ways to calculate small but strong deviations beyond using simple moving averages.
Some options I've explored include EMA, z-scores, step detection, and anomaly detection. I would greatly appreciate any insights or suggestions.


text:
rna seq matrix design hybrids
I am struggling to get my head around the correct design to use for my experiment.
I have RNAseq of 2 pure species (s1, s2) and of the hybrid (h). I aligned all three genotypes to the s1 genomes. I am interested in the comparison of hybrids vs s1 (maternal genome). But I want to remove a potential bias from reads from s2 chromosome being harder to align to s1. (The hybrid genome is composed of both s1 and s2.) I would like to correct my DE test by removing potential effect s2 vs s1. (I aligned s2 to s1 and got the count.)
But my issue is, I fail to have a full rank design.
I do not have any strong feelings about the library (deseq2, edgeR...).
Following the tutorial
below is the condition matrix, as you can see if I want to make the contrast ~ s2 + condition. The matrix is not full rank.
I believe I could write my design as ~s2 + h - s1 ?? But I don't know how to use deseq2 for this approach.
Any approach/tips to achieve my goal here? Is it event possible? what

sample
condition
s2

1
s1
0

2
s1
0

3
s1
0

4
s2
1

5
s2
1

6
s2
1

7
h
0

8
h
0

9
h
0

any help is warmly appreciated!

I think you are over thinking the problem. You need to classify the genes of interest and to do that a classifier such as eggNOG is used. This will remove any bias. I understand eggNOG because it's based on protein phylogenetics and phylogenetics is what I do and I can assure you its very solid and the algorithm nicely constructed.
DESeq2 (I ain't a huge fan BTW, but avoid controversy) will happily accept eggNOG classification, here
Then its simply three tests

S1 vs S2
S1 vs hybrid
S2 vs hybrid

I personally would perform unsupervised learning, particularly via Scanpy first on S1, S2 and hybrid separately. This will give a rough indication of the key difference of expression profiles between S1 and S2 and see where the hydrid occurs.


text:
Generating synthetic RNA-Seq data for performance assessment of new algorithm
I want to generate synthetic miRNA data for the performance assessment of the newly developed algorithm designed for miRNA data analysis. There are several synthetic sequence simulators available in the published literature, such as wgsim, dwgsim, ART, inSillicoseq, mason,  NEAT, BEER, polyester, RSEM, SimSeq, etc. All of these tools are either for DNA or transcriptome profiles. Can anybody please suggest the following:

Can synthetic DNA sequence simulators (ART, wgsim, dwgsim) generate synthetic miRNA data?
Some tools, such as CAMPAREE, RSEM, and RNA-Seq simulator, generate RNA-seq data, but I think these tools generate mRNA data, and I am unsure whether I can use them for my miRNA data analysis.

Please suggest. I very much appreciate any help in this regard. Thanks


text:
How do I do "protein culling"?
I am working on a novel protein secondary structure detection algorithm using machine learning.
As far as I understand, online protein databanks (e.g., RCSB, etc.) have a huge collection of proteins (several millions). All proteins are not useful for research purposes as they contain redundant, duplicate, or incomplete information.
So, we have to choose the minimum number of proteins so that we have enough proteins to cover all the required information for our specific purpose.
This is called "protein culling."
Google search shows that there were some servers like PISCES, ProDaMa, etc. However, as usual, all of them are down.
Now the question is, how do I do protein-culling myself in the context of protein secondary structure detection?
N.B. Proteins with similar folds or secondary structure elements can be considered redundant, even if their sequence identity is low.

How good are your data engineering skills? Removing duplicate, redundant or incomplete stuff isn't hard and a database of several million is honestly pretty small. You can probably do this in pandas2 definitely polars, but after that it gets more specialised.

@M__, I work on secondary and tertiary structures. I need a minimum set of individual PDB files.

@M__, PIESCES server had 5022 PDB files.

It's still applicable to data engineering for the criteria you are considering. I agree that it may not fit inside pandas2. If you were considering genetic distance ... thats more complicated.

This depends entirely on what, exactly, you are trying to do. Please [edit] your question and specify the biological question you are trying to answer.

@terdon, I am working on a novel protein secondary and tertiary structure detection algorithm using machine learning.

OK, please add that to your question but, again, what do you consider not useful here? Presumably, you don't care so much about sequence identity as about tertiary structure similarity, right? So what would you like to cull and what would you like to keep?

@terdon, Proteins with similar folds or secondary or tertiary structure elements can be considered redundant, even if their sequence identity is low.

Have you considered looking at this in terms of protein domains instead of proteins? Say starting with the Pfam data and getting representative examples of each domain? (random shot in the dark here, I have no experience with tertiary structure analysis)

@terdon, I am a beginner at level zero in Bioinformatics.

If you were willing to download the data yourself, you could use a tool like FoldSeek to cluster on various thresholds to remove redundancy. Though I realise this is tertiary structure, rather than secondary structure, but it still could work with some tweaking.


text:
How to find a database for proteins that bind to different molecules and have different binding sites?
I'm doing a study in which I'm trying to analyze the effect of molecules (ligands) on binding sites of proteins.
Therefore, I'm trying to find a database in which the proteins bind to multiple different molecules, and as a result have different binding sites. The database that I'm working with currently is from UniProt, and the problem with it is that it has proteins that bind to only a small class of molecules, and hence the binding site is similar for each protein-ligand pair.
For reference, here is the database that I'm using right now: https://bit.ly/ProteinLigandBindingSite
Can someone help me find a database which contains proteins that bind to different molecules, and hence have different binding sites?
Thank you!

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.


text:
How to perform meaningful Gene set erichment analysis or otherwise find broader themes/functions in different cell populations?
(For context: I'm somewhat new to bioinformatic analyses but am mostly comfortable with R)
I have identified transcriptomically and morphologically different cancer cell populations within a patient (scRNAseq, spatial transcriptomics, IHC). They are quite distinct, with hundreds of genes as differentially expressed genes (DEGs) even after high tresholds. It's quite interesting to describe, but I am looking to find better ways to represent their features. Of course I am attaching DEG tables as supplementary data etc. but I thought that gene set enrichment analysis could give broad terms that would a) help me understand which DEGs might contribute together to biological processes and b) describe the bigger picture.
So, I tried running a few terms (GO:BP, GO:CC, GO:MF, ...) from the GSEA msigdb website on the DEGs from scRNAseq, but have come out a bit frustrated. A lot of terms make sense and reflect the changes, but have large overlap. For example, the 3 top GO:BP terms (Extracellular matrix organization", "Extracellular structure organization", "External encapsulating structure organziation" are scored identically, with identical DEGs as the genes that are part of their gene sets. This is repeated further down the list with subsets of this gene set. Now, reporting these doesn't actually do what I want, which is give a global overview (but the information is heavily redundant). I could just pick one of the 3 redundant ones and skip the others, but then I am worried my presentation would be biased.
Any advice on how to properly do gene set enrichment analysis to find broader terms to describe and understand what's going on? Or alternative ways that lead to the goal of broader description and understanding? I would be very grateful.
Thank you!

Additionally there are plenty of tools to summarize long lists of GO terms by removing redundancy (which seems to be your concern), such as REVIGO (web) or rrvgo (R+Bioconductor). Disclaimer: I'm the author of the latter.
Such tools objectively pick the best representative term of each set of redundant terms, eg. by considering the GO term size (in genes), observed/expected ratio, p-value, etc.

It is totally acceptable in your main figure to display only a selected part of your GO term, as long as the reader is aware of this action (e.g., by including a legend specifying "selected pathways enriched"). Additionally, ensure that either readers have sufficient information to replicate your analyses or provide the full result table as supplementary information. By adhering to these points, you are transparent on what you did, allowing others to replicate the analyses independently. Purposely excluding top hit pathways that contradict your theory could be considered misleading.
If you are looking for broader terms, have you tried to a use slim GO? Those are GO with only a subset of the term. You could also look at the GO graph for your term and select the one higher in the graph.
There are also plenty of tools out there to display GO results-graph. GOplot, GOnet are the first hit I got on google.