text:
How can I resolve "missing name of population" message in PGDSpider?
I want to use Arlequin for AMOVA. I have my FASTA files which I want to convert into .arp files which are recognizable by arlequin. When I try to convert it using Population Genetics Data (format) Spider, or PGDSpider by Lisher and Excoffier (2018).
The error message says the name of my populations is missing.
Can anyone tell me what's wrong with my FASTA files? Check the following example:
>CcHr1
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT
ACTAATTCGTCTTGAGTTAGGTACTTCTGGTGTCTTGACTGATGAGCATTTTTATAATGTAATTGTAACA
GCCCACGCTTTTGTAATGATTTTTTTTATGGTTATACCGATTATAATTGGTGGTTTTGGAAATTGAATAG
TTCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGGATAAATAATATAAGATTTTGATTATTACC
TCCTGCTTTTGTACTACTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGGACTGGCTGAACTGTTTAT
CCTCCACTTAGGAGGGTAACAGGACATAGTGGGGCTTCTGTAGACTTGGCTATTTTTTCCTTACATCTCG
CGGGCATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAGCGGATAAGGCTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTACTATCACTT
CCAGTTCTTGCAGGGGCAATTACGATGCTCTTAACTGATCGTAATTTTAGTACTTCATTTTTTGATCCAG
>CcMath1
ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT
ATTAATTCGTCTTGAGTTAGGTACTTCTGGTATCTTGACTGATGAGCATTTTTATAACGTAATTGTAACA
GCCCACGCTTTTGTAATAATTTTTTTTATGGTTATGCCAATTATAATTGGTGGTTTTGGAAATTGAATAG
TCCCTTTATTAATTGGTGCTCCAGACATAAGGTTCCCTCGAATAAATAATATAAGATTTTGACTATTACC
TCCTGCGTTTGTACTGCTGATTAGTAGAAGCCTTGTTGAGGGAGGGGCTGGAACTGGCTGAACTGTTTAT
CCTCCACTTAGAAGGGTAACAGGGCATAGTGGGGCTTCTGTAGACTTGGCGATCTTTTCCTTGCATCTCG
CGGGTATATCATCAATTTTGGGGGCAATTAATTTTATTACAACGATCTTCAATATGCGGTCCCCTGGGTT
AACACTCGAACGGATAAGGTTATTTGTATGGTCTATTTTAGTAACAGTTTTTTTATTGTTATTGTCACTT
CCAGTTCTTGCAGGGGCTATTACGATACTCTTAACTGATCGTAATTTTAGTACTTCGTTTTTTGACCCAG

As you can see I have to populations but it only names it as pop1.

Is this something I can do inside arlequin? Because so far I have only though of manipulating the flat file into distinguishing populations. This is a very tedious process I would like to avoid

You could try specifying the population keyword in the record header. From the PGDSpider manual (version 2.1.1.5 (May 2018)), at the bottom of page 67:

Population definition in header line: String after “population:” until
the first space is parsed as population name e.g.: >sample1 population:pop1 gene:HMA4

A good response.

It's important to understand what Arlequin is doing for MANOVA. It will compare two or more populations to assess whether there is a difference in the allele frequency between the populations. It will also assess "interaction" (interpreting that in "gene flow" needs some thinking).
The basic theory is that if two or more populations cannot be distinguished under a formal parametric test there must be very strong gene flow supporting the lack of population structure between these populations.
Geographic foci Population designation is usually done between geographic foci, thus population A is usually geographic location (focus) A and population B is a geographically separated region (focus). Thus the fasta header needs to know the geographic region of the sample. Thats it.
Population info The issue is that information is required within the fasta header. Thus population delineation is therefore essential. This is confirmed in the PGDSpider manual below the line for Arlequin: which you don't need to do because thats what the conversion does.
Arlequin therefore needs to group each member of the population into a (usually geographically seperated) population to perform MANOVA. It will require this information specifying on the fasta header line. At the moment the data is saying there is a single population, it thinks all the sequences are from the same geographic location.
Procedure Thus it's likely simply a case of assigning two or more broad geographic foci for your samples. Then going through the fasta file and designating each sequence into one of the geographic focus under the keyword "population".
Just modifying @Steve's reading of the manual
> CcHr1 population:geoFocusA gene:X
ATAAAGATATTGGTACGTTGTATATATTATTGGGTGTGTGGTGTGGTATAGTTGGTACAGGATTATCATT

>CcMath1 population:geoFocusB gene:X
 ATAAAGATATTGGTACGTTATATATATTATTGGGTGTGTGGTGTGGAATAGTTGGCACAGGGTTGTCATT

Then run PGDSpider
java -Xmx1024m -Xms512m -jar PGDSpider2-cli.jar -inputfile examples\example_Structure.fasta -inputformat FASTA -outputfile examples\output_Arlequin.arp -outputformat ARLEQUIN -spid examples\Structure_Arlequin.spid

You did attain the critical anomaly in the sequence format,

each line of a sequence should have fewer than 80 characters

Your rows have 70 residues.

Genetic structure (only required for AMOVA):

The genetic structure specifies the hierarchical genetic structure of the samples. It is possible to define groups of populations.
o start of the subsection: [[Structure]]
o name for the genetic structure (string within ” ”):  StructureName = “name”
o number of groups defined in the structure (int value): NbGroups = 5
o group definitions (list containing the names of the samples belonging to the group,
entered within braces “{ }”):
Id1 2 ACTCGGGTTCGCGCGC # the first pseudo-haplotype ACTCGGGCTCACGCGC # the second pseudo-haplotype

NbGroups=2
Group ={
  population1
}
Group ={
  population2
  population3
}

Advantage The advantage of MANOVA is you can have LOADS of populations and it will assess which have high levels of gene flow between them. Statistically it is fantastic for loads of different reasons (quite complicated reasons). The disadvantage of output of MANOVA is there will likely be  significant differences between many of the populations - possibly all of them. This is a result of the power of parametric statistics, but may not be particularly helpful for the analysis. If that is the outcome switching to Fst is the next part of the analysis (and a different question). Goodluck and you can always post back if you get stuck.


text:
How can I fix "Ill-defined genetic structure" Error message in Arlequin?
I am performing an  AMOVA analysis via Arlequin.
At this point, I have:

collected sequences of mtDNA from the cox1 gene,
aligned the different genes using clustalX
converted the files into .arp format.

I then altered the structure of my populations in Arlequin so that all the different populations are in one group. When I tried to run AMOVA I get the error message

"Ill-Defined genetic structure".

Any clues as to what I have to adjust?

If one group is one population this will not work with Arlequin AMOVA. There must be more than one population (I might have misread the question however).
There could several reasons for this output. The most likely is there are insufficient SNPs within the cox1 locus for MANOVA to identify a population structure with sufficient statistical robustness. What it requires are a minimum of 3 independent alleles (SNPs), really 3 genotypes. If the SNPs are linked on the same bit of DNA that doesn't count - thats just one genotype. What I suspect is happening is there is a dearth (not enough) genetic information within cox1 to perform MANOVA.
There could be other explanations, but this is the first issue to address. This isn't just Arlequin it is standard theory for any allele frequency based calculations and why this is important is more complicated to explain.

Please note, it is the run-up to Christmas .. okay put Christmas on hold. To respond to the comments ...
Firstly, if you hold the data set it is very easy to identify the number of genotypes. There are past posts on this herein. I am not inclined to describe alternative approaches to an algorithm I've still not yet published :-) ... however it is really easy and needs an empirical test. It is good to assess number of genotypes per se Arlequin will do this analysis ... finding that within the output I dunno.
Extending loci is always good and for mtDNA I would not recommend with 16S nor cox2. rRNA should be avoided in allele frequency calculations, its okay for taxonomic assignment (bacteria rather than mtDNA). cox2 is a small locus. I would recommend for mtDNA

cobB, formerly called cytB
intergenic region of mtDNA

cobB is a good target, it has a lot more selection on it for within species analysis than cox1. The intergenic region is rarely used, its a good locus because of the higher rates of mutation therein.

Would you suggest using a different genetic marker for this type of analyses? I do not have access to SNPs so the next best thing would be other mtDNA genes like COX2 and rRNA.

@Nickmofoe answered ... Merry Christmas BTW.

Alright so: The COX1 locus ended up not being the issue. I made some corrections to the flat files and they seemed to run with no issues. I got the results and everything seems to check out. I concluded that the issue was not the locus itself but rather the incredibly complicated .arp format. Thanks for the advice though! Much appreciated and merry Christmas!!!

Merry Christmas to you @nickmofoe! The  "complicated format" was why Arlequin was never popular, because everybody felt it was unnecessarily complicated. This why your original question on PGDSpider was interesting because it as overdue requirement to make Arlequin accessible.

source: Stack Exchange

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

sort: votes


text:
Error installing VEP on mac
I have tried multiple different ways to get Ensemble's Variant Effect Predictor (VEP) to properly install on my mac, but I keep getting same error.
In /Users/username/ensembl-vep/htslib
make: Nothing to be done for `all'.
 - unpacking ./Bio/tmp/biodbhts.zip to ./Bio/tmp/
./Bio/tmp/Bio-DB-HTS-2.11 - moving files to ./biodbhts
 - making Bio::DB:HTS
clang-4.0: error: invalid version number in '-mmacosx-version-min=12.4'
Warning: ExtUtils::CBuilder not installed or no compiler detected
Proceeding with configuration, but compilation may fail during Build

Created MYMETA.yml and MYMETA.json
Creating new 'Build' script for 'Bio-DB-HTS' version '2.11'
Building Bio-DB-HTS
Error: no compiler detected to compile 'lib/Bio/DB/HTS.c'.  Aborting
ERROR: Shared Bio::DB:HTS library not found

--> When I check for the version of clang, I always get version 12, not sure where clang-4.0 is. I think error is specifically not correct path to current version of C++.

Do you have xcode utilities installed on your mac? What else have you tried?

You can try and get Perl's ExUtils::Builder, described here correctly installed and working, for example
perl -MCPAN -e shell
force install Bundle::CPAN

You might want to use brew to make the build safe.
I don't do this stuff, I use OSX, I would install ensemble's VEP via conda
conda install -c bioconda ensembl-vep

You might want to (good idea) create a specific environment and obviously download Anaconda or miniconda.
conda create --name vep
conda activate vep

It will work, the output is below. I didn't proceed with the installation because of the shear amount of stuff downloaded into a given conda environment.  There's everything in there but the kitchen sink, like there is mafft, t-coffee, clustalw, muscle and clustalo required in the download. I am honestly clueless because they are all alignment algorithms. Anyway I honestly wouldn't chase the MCPAN stuff, in this instance I definitely recommend conda. Finally, you honestly you don't want all this extra stuff on your hard drive when you've finished just delete the entire environment.
The final point I'd like to make is muscle5 (which is also in there) has recently significantly been modified making all other aligners obsolete on the runtime issue, but has to be run under the super5 mode.
## Package Plan ##

  environment location: /opt/anaconda3/envs/XXX # environment was edited  out

  added / updated specs:
    - ensembl-vep

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    argtable2-2.13             |    h1de35cc_1001         2.7 MB  conda-forge
    clustalo-1.2.4             |       h9722bc1_5         331 KB  bioconda
    clustalw-2.1               |       hcd10b59_7         326 KB  bioconda
    ensembl-vep-108.2          | pl5321h4a94de4_0         3.1 MB  bioconda
    expat-2.4.9                |       he9d5cce_0         127 KB
    htslib-1.14                |       hc38c3fb_0         2.0 MB  bioconda
    libdeflate-1.7             |       h9ed2024_5          52 KB
    libgfortran-5.0.0          |11_3_0_hecd8cb5_28         142 KB
    libgfortran5-11.3.0        |      h9dfd629_28         1.4 MB
    libzlib-1.2.12             |       hfd90126_3          65 KB  conda-forge
    mafft-7.508                |       ha5712d3_0         8.7 MB  bioconda
    muscle-5.1                 |       hb339e23_1         298 KB  bioconda
    mysql-connector-c-6.1.11   |       hccea1a4_1         1.2 MB
    paml-4.9                   |       hb4d813b_6         1.2 MB  bioconda
    perl-algorithm-diff-1.1903 | pl5321hdfd78af_3          36 KB  bioconda
    perl-base-2.23             | pl5321hdfd78af_2          12 KB  bioconda
    perl-bio-asn1-entrezgene-1.73| pl5321hdfd78af_3          24 KB  bioconda
    perl-bio-coordinate-1.007001| pl5321hdfd78af_3          36 KB  bioconda
    perl-bio-db-hts-3.01       | pl5321hc7384eb_4         116 KB  bioconda
    perl-bio-featureio-1.6.905 | pl5321hdfd78af_4          41 KB  bioconda
    perl-bio-samtools-1.43     | pl5321h1f540d2_3         161 KB  bioconda
    perl-bio-tools-phylo-paml-1.7.3| pl5321hdfd78af_3          58 KB  bioconda
    perl-bio-tools-run-alignment-clustalw-1.7.4| pl5321hdfd78af_3          25 KB  bioconda
    perl-bio-tools-run-alignment-tcoffee-1.7.4| pl5321hdfd78af_5          25 KB  bioconda
    perl-bioperl-1.7.8         |       hdfd78af_0           8 KB  bioconda
    perl-bioperl-core-1.7.8    | pl5321h527b516_0         1.8 MB  bioconda
    perl-bioperl-run-1.007003  | pl5321hdfd78af_0         340 KB  bioconda
    perl-business-isbn-3.007   | pl5321hdfd78af_0          18 KB  bioconda
    perl-business-isbn-data-20210112.006| pl5321hdfd78af_0          21 KB  bioconda
    perl-capture-tiny-0.48     | pl5321hdfd78af_2          19 KB  bioconda
    perl-class-data-inheritable-0.09| pl5321hdfd78af_0          10 KB  bioconda
    perl-constant-1.33         | pl5321hdfd78af_2          12 KB  bioconda
    perl-data-dumper-2.183     | pl5321ha5712d3_1          38 KB  bioconda
    perl-dbd-mysql-4.046       | pl5321hcd10b59_4          94 KB  bioconda
    perl-dbi-1.643             | pl5321ha5712d3_1         602 KB  bioconda
    perl-devel-stacktrace-2.04 | pl5321hdfd78af_1          17 KB  bioconda
    perl-digest-hmac-1.04      | pl5321hdfd78af_0          10 KB  bioconda
    perl-digest-md5-2.58       | pl5321ha5712d3_1          18 KB  bioconda
    perl-encode-locale-1.05    | pl5321hdfd78af_7          13 KB  bioconda
    perl-exception-class-1.45  | pl5321hdfd78af_0          21 KB  bioconda
    perl-file-listing-6.15     | pl5321hdfd78af_0          14 KB  bioconda
    perl-file-slurp-tiny-0.004 | pl5321hdfd78af_2          10 KB  bioconda
    perl-file-sort-1.01        | pl5321hdfd78af_3          20 KB  bioconda
    perl-file-spec-3.48_01     | pl5321hdfd78af_2           6 KB  bioconda
    perl-getopt-long-2.54      | pl5321hdfd78af_0          31 KB  bioconda
    perl-html-parser-3.80      | pl5321hb3ea9af_0          67 KB  bioconda
    perl-html-tagset-3.20      | pl5321hdfd78af_4          14 KB  bioconda
    perl-http-cookies-6.10     | pl5321hdfd78af_0          22 KB  bioconda
    perl-http-daemon-6.14      | pl5321hdfd78af_0          20 KB  bioconda
    perl-http-date-6.05        | pl5321hdfd78af_0          13 KB  bioconda
    perl-http-message-6.36     | pl5321hdfd78af_0          54 KB  bioconda
    perl-http-negotiate-6.01   | pl5321hdfd78af_4          14 KB  bioconda
    perl-io-html-1.004         | pl5321hdfd78af_0          16 KB  bioconda
    perl-io-socket-ssl-2.074   | pl5321hdfd78af_0         162 KB  bioconda
    perl-io-string-1.08        | pl5321hdfd78af_4          12 KB  bioconda
    perl-io-tty-1.16           | pl5321ha5712d3_1          30 KB  bioconda
    perl-ipc-run-20200505.0    | pl5321hdfd78af_0          76 KB  bioconda
    perl-libwww-perl-6.67      | pl5321hdfd78af_0          97 KB  bioconda
    perl-libxml-perl-0.08      | pl5321hdfd78af_3          40 KB  bioconda
    perl-lwp-mediatypes-6.04   | pl5321hdfd78af_1          24 KB  bioconda
    perl-mime-base64-3.16      | pl5321ha5712d3_2          16 KB  bioconda
    perl-net-http-6.22         | pl5321hdfd78af_0          21 KB  bioconda
    perl-net-ssleay-1.92       | pl5321h5d619eb_1         299 KB  bioconda
    perl-ntlm-1.09             | pl5321hdfd78af_5          16 KB  bioconda
    perl-perlio-gzip-0.20      | pl5321h1f540d2_3          16 KB  bioconda
    perl-sereal-4.019          | pl5321hdfd78af_0          11 KB  bioconda
    perl-sereal-decoder-4.025  | pl5321ha5712d3_1         156 KB  bioconda
    perl-sereal-encoder-4.025  | pl5321ha5712d3_1         365 KB  bioconda
    perl-set-intervaltree-0.12 | pl5321h9722bc1_2          31 KB  bioconda
    perl-socket-2.027          | pl5321ha5712d3_3          29 KB  bioconda
    perl-sub-uplevel-0.2800    | pl5321ha5712d3_4          14 KB  bioconda
    perl-test-deep-1.130       | pl5321hdfd78af_0          37 KB  bioconda
    perl-test-differences-0.69 | pl5321hdfd78af_0          19 KB  bioconda
    perl-test-exception-0.43   | pl5321hdfd78af_3          15 KB  bioconda
    perl-test-harness-3.44     | pl5321hdfd78af_0          85 KB  bioconda
    perl-test-most-0.38        | pl5321hdfd78af_0          20 KB  bioconda
    perl-test-warn-0.36        | pl5321hdfd78af_2          14 KB  bioconda
    perl-text-csv-2.01         | pl5321hdfd78af_0          70 KB  bioconda
    perl-text-diff-1.45        | pl5321hdfd78af_1          23 KB  bioconda
    perl-time-local-1.30       | pl5321hdfd78af_0          12 KB  bioconda
    perl-timedate-2.33         | pl5321hdfd78af_2          31 KB  bioconda
    perl-tree-dag_node-1.32    | pl5321hdfd78af_0          43 KB  bioconda
    perl-try-tiny-0.31         | pl5321hdfd78af_1          19 KB  bioconda
    perl-uri-5.12              | pl5321hdfd78af_0          59 KB  bioconda
    perl-url-encode-0.03       | pl5321h527b516_0          12 KB  bioconda
    perl-www-robotrules-6.02   | pl5321hdfd78af_4          14 KB  bioconda
    perl-xml-dom-1.46          | pl5321hdfd78af_1          59 KB  bioconda
    perl-xml-dom-xpath-0.14    | pl5321hdfd78af_2          12 KB  bioconda
    perl-xml-parser-2.44       | pl5321h50c80e2_8         162 KB  bioconda
    perl-xml-regexp-0.04       | pl5321hdfd78af_3          10 KB  bioconda
    perl-xml-xpathengine-0.14  | pl5321hdfd78af_3          28 KB  bioconda
    t-coffee-11.00.8cbe486     |       h26a2512_0         1.1 MB  bioconda
    unzip-6.0                  |       h3b3adb9_0          71 KB
    ------------------------------------------------------------
                                           Total:        28.7 MB

The following NEW packages will be INSTALLED:

  argtable2          conda-forge/osx-64::argtable2-2.13-h1de35cc_1001
  blast              bioconda/osx-64::blast-2.12.0-h0370960_3
  clustalo           bioconda/osx-64::clustalo-1.2.4-h9722bc1_5
  clustalw           bioconda/osx-64::clustalw-2.1-hcd10b59_7
  ensembl-vep        bioconda/noarch::ensembl-vep-108.2-pl5321h4a94de4_0
  entrez-direct      bioconda/osx-64::entrez-direct-16.2-h193322a_1
  expat              pkgs/main/osx-64::expat-2.4.9-he9d5cce_0
  htslib             bioconda/osx-64::htslib-1.14-hc38c3fb_0
  libdeflate         pkgs/main/osx-64::libdeflate-1.7-h9ed2024_5
  libgfortran        pkgs/main/osx-64::libgfortran-5.0.0-11_3_0_hecd8cb5_28
  libgfortran5       pkgs/main/osx-64::libgfortran5-11.3.0-h9dfd629_28
  libidn2            pkgs/main/osx-64::libidn2-2.3.2-h9ed2024_0
  libunistring       pkgs/main/osx-64::libunistring-0.9.10-h9ed2024_0
  libzlib            conda-forge/osx-64::libzlib-1.2.12-hfd90126_3
  llvm-openmp        pkgs/main/osx-64::llvm-openmp-14.0.6-h0dcd299_0
  mafft              bioconda/osx-64::mafft-7.508-ha5712d3_0
  muscle             bioconda/osx-64::muscle-5.1-hb339e23_1
  mysql-connector-c  pkgs/main/osx-64::mysql-connector-c-6.1.11-hccea1a4_1
  paml               bioconda/osx-64::paml-4.9-hb4d813b_6
  pcre               pkgs/main/osx-64::pcre-8.45-h23ab428_0
  perl               pkgs/main/osx-64::perl-5.32.1-0_h435f0c2_perl5
  perl-algorithm-di~ bioconda/noarch::perl-algorithm-diff-1.1903-pl5321hdfd78af_3
  perl-archive-tar   bioconda/noarch::perl-archive-tar-2.40-pl5321hdfd78af_0
  perl-base          bioconda/noarch::perl-base-2.23-pl5321hdfd78af_2
  perl-bio-asn1-ent~ bioconda/noarch::perl-bio-asn1-entrezgene-1.73-pl5321hdfd78af_3
  perl-bio-coordina~ bioconda/noarch::perl-bio-coordinate-1.007001-pl5321hdfd78af_3
  perl-bio-db-hts    bioconda/osx-64::perl-bio-db-hts-3.01-pl5321hc7384eb_4
  perl-bio-featureio bioconda/noarch::perl-bio-featureio-1.6.905-pl5321hdfd78af_4
  perl-bio-samtools  bioconda/osx-64::perl-bio-samtools-1.43-pl5321h1f540d2_3
  perl-bio-tools-ph~ bioconda/noarch::perl-bio-tools-phylo-paml-1.7.3-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-clustalw-1.7.4-pl5321hdfd78af_3
  perl-bio-tools-ru~ bioconda/noarch::perl-bio-tools-run-alignment-tcoffee-1.7.4-pl5321hdfd78af_5
  perl-bioperl       bioconda/noarch::perl-bioperl-1.7.8-hdfd78af_0
  perl-bioperl-core  bioconda/osx-64::perl-bioperl-core-1.7.8-pl5321h527b516_0
  perl-bioperl-run   bioconda/noarch::perl-bioperl-run-1.007003-pl5321hdfd78af_0
  perl-business-isbn bioconda/noarch::perl-business-isbn-3.007-pl5321hdfd78af_0
  perl-business-isb~ bioconda/noarch::perl-business-isbn-data-20210112.006-pl5321hdfd78af_0
  perl-capture-tiny  bioconda/noarch::perl-capture-tiny-0.48-pl5321hdfd78af_2
  perl-carp          bioconda/noarch::perl-carp-1.38-pl5321hdfd78af_4
  perl-class-data-i~ bioconda/noarch::perl-class-data-inheritable-0.09-pl5321hdfd78af_0
  perl-common-sense  bioconda/noarch::perl-common-sense-3.75-pl5321hdfd78af_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-bzip2-2.201-pl5321h9722bc1_0
  perl-compress-raw~ bioconda/osx-64::perl-compress-raw-zlib-2.105-pl5321h9722bc1_0
  perl-constant      bioconda/noarch::perl-constant-1.33-pl5321hdfd78af_2
  perl-data-dumper   bioconda/osx-64::perl-data-dumper-2.183-pl5321ha5712d3_1
  perl-dbd-mysql     bioconda/osx-64::perl-dbd-mysql-4.046-pl5321hcd10b59_4
  perl-dbi           bioconda/osx-64::perl-dbi-1.643-pl5321ha5712d3_1
  perl-devel-stackt~ bioconda/noarch::perl-devel-stacktrace-2.04-pl5321hdfd78af_1
  perl-digest-hmac   bioconda/noarch::perl-digest-hmac-1.04-pl5321hdfd78af_0
  perl-digest-md5    bioconda/osx-64::perl-digest-md5-2.58-pl5321ha5712d3_1
  perl-encode        bioconda/osx-64::perl-encode-3.19-pl5321ha5712d3_1
  perl-encode-locale bioconda/noarch::perl-encode-locale-1.05-pl5321hdfd78af_7
  perl-exception-cl~ bioconda/noarch::perl-exception-class-1.45-pl5321hdfd78af_0
  perl-exporter      bioconda/noarch::perl-exporter-5.72-pl5321hdfd78af_2
  perl-exporter-tiny bioconda/noarch::perl-exporter-tiny-1.002002-pl5321hdfd78af_0
  perl-extutils-mak~ conda-forge/noarch::perl-extutils-makemaker-7.64-pl5321hd8ed1ab_0
  perl-file-listing  bioconda/noarch::perl-file-listing-6.15-pl5321hdfd78af_0
  perl-file-slurp-t~ bioconda/noarch::perl-file-slurp-tiny-0.004-pl5321hdfd78af_2
  perl-file-sort     bioconda/noarch::perl-file-sort-1.01-pl5321hdfd78af_3
  perl-file-spec     bioconda/noarch::perl-file-spec-3.48_01-pl5321hdfd78af_2
  perl-getopt-long   bioconda/noarch::perl-getopt-long-2.54-pl5321hdfd78af_0
  perl-html-parser   bioconda/osx-64::perl-html-parser-3.80-pl5321hb3ea9af_0
  perl-html-tagset   bioconda/noarch::perl-html-tagset-3.20-pl5321hdfd78af_4
  perl-http-cookies  bioconda/noarch::perl-http-cookies-6.10-pl5321hdfd78af_0
  perl-http-daemon   bioconda/noarch::perl-http-daemon-6.14-pl5321hdfd78af_0
  perl-http-date     bioconda/noarch::perl-http-date-6.05-pl5321hdfd78af_0
  perl-http-message  bioconda/noarch::perl-http-message-6.36-pl5321hdfd78af_0
  perl-http-negotia~ bioconda/noarch::perl-http-negotiate-6.01-pl5321hdfd78af_4
  perl-io-compress   bioconda/osx-64::perl-io-compress-2.201-pl5321h9722bc1_0
  perl-io-html       bioconda/noarch::perl-io-html-1.004-pl5321hdfd78af_0
  perl-io-socket-ssl bioconda/noarch::perl-io-socket-ssl-2.074-pl5321hdfd78af_0
  perl-io-string     bioconda/noarch::perl-io-string-1.08-pl5321hdfd78af_4
  perl-io-tty        bioconda/osx-64::perl-io-tty-1.16-pl5321ha5712d3_1
  perl-io-zlib       bioconda/noarch::perl-io-zlib-1.11-pl5321hdfd78af_0
  perl-ipc-run       bioconda/noarch::perl-ipc-run-20200505.0-pl5321hdfd78af_0
  perl-json          bioconda/noarch::perl-json-4.10-pl5321hdfd78af_0
  perl-json-xs       bioconda/osx-64::perl-json-xs-2.34-pl5321hcd10b59_5
  perl-libwww-perl   bioconda/noarch::perl-libwww-perl-6.67-pl5321hdfd78af_0
  perl-libxml-perl   bioconda/noarch::perl-libxml-perl-0.08-pl5321hdfd78af_3
  perl-list-moreuti~ bioconda/noarch::perl-list-moreutils-0.430-pl5321hdfd78af_0
  perl-list-moreuti~ bioconda/osx-64::perl-list-moreutils-xs-0.430-pl5321ha5712d3_1
  perl-lwp-mediatyp~ bioconda/noarch::perl-lwp-mediatypes-6.04-pl5321hdfd78af_1
  perl-mime-base64   bioconda/osx-64::perl-mime-base64-3.16-pl5321ha5712d3_2
  perl-net-http      bioconda/noarch::perl-net-http-6.22-pl5321hdfd78af_0
  perl-net-ssleay    bioconda/osx-64::perl-net-ssleay-1.92-pl5321h5d619eb_1
  perl-ntlm          bioconda/noarch::perl-ntlm-1.09-pl5321hdfd78af_5
  perl-parent        bioconda/noarch::perl-parent-0.236-pl5321hdfd78af_2
  perl-pathtools     bioconda/osx-64::perl-pathtools-3.75-pl5321ha5712d3_3
  perl-perlio-gzip   bioconda/osx-64::perl-perlio-gzip-0.20-pl5321h1f540d2_3
  perl-scalar-list-~ bioconda/osx-64::perl-scalar-list-utils-1.62-pl5321ha5712d3_0
  perl-sereal        bioconda/noarch::perl-sereal-4.019-pl5321hdfd78af_0
  perl-sereal-decod~ bioconda/osx-64::perl-sereal-decoder-4.025-pl5321ha5712d3_1
  perl-sereal-encod~ bioconda/osx-64::perl-sereal-encoder-4.025-pl5321ha5712d3_1
  perl-set-interval~ bioconda/osx-64::perl-set-intervaltree-0.12-pl5321h9722bc1_2
  perl-socket        bioconda/osx-64::perl-socket-2.027-pl5321ha5712d3_3
  perl-sub-uplevel   bioconda/osx-64::perl-sub-uplevel-0.2800-pl5321ha5712d3_4
  perl-test-deep     bioconda/noarch::perl-test-deep-1.130-pl5321hdfd78af_0
  perl-test-differe~ bioconda/noarch::perl-test-differences-0.69-pl5321hdfd78af_0
  perl-test-excepti~ bioconda/noarch::perl-test-exception-0.43-pl5321hdfd78af_3
  perl-test-harness  bioconda/noarch::perl-test-harness-3.44-pl5321hdfd78af_0
  perl-test-most     bioconda/noarch::perl-test-most-0.38-pl5321hdfd78af_0
  perl-test-warn     bioconda/noarch::perl-test-warn-0.36-pl5321hdfd78af_2
  perl-text-csv      bioconda/noarch::perl-text-csv-2.01-pl5321hdfd78af_0
  perl-text-diff     bioconda/noarch::perl-text-diff-1.45-pl5321hdfd78af_1
  perl-time-local    bioconda/noarch::perl-time-local-1.30-pl5321hdfd78af_0
  perl-timedate      bioconda/noarch::perl-timedate-2.33-pl5321hdfd78af_2
  perl-tree-dag_node bioconda/noarch::perl-tree-dag_node-1.32-pl5321hdfd78af_0
  perl-try-tiny      bioconda/noarch::perl-try-tiny-0.31-pl5321hdfd78af_1
  perl-types-serial~ bioconda/noarch::perl-types-serialiser-1.01-pl5321hdfd78af_0
  perl-uri           bioconda/noarch::perl-uri-5.12-pl5321hdfd78af_0
  perl-url-encode    bioconda/osx-64::perl-url-encode-0.03-pl5321h527b516_0
  perl-www-robotrul~ bioconda/noarch::perl-www-robotrules-6.02-pl5321hdfd78af_4
  perl-xml-dom       bioconda/noarch::perl-xml-dom-1.46-pl5321hdfd78af_1
  perl-xml-dom-xpath bioconda/noarch::perl-xml-dom-xpath-0.14-pl5321hdfd78af_2
  perl-xml-parser    bioconda/osx-64::perl-xml-parser-2.44-pl5321h50c80e2_8
  perl-xml-regexp    bioconda/noarch::perl-xml-regexp-0.04-pl5321hdfd78af_3
  perl-xml-xpatheng~ bioconda/noarch::perl-xml-xpathengine-0.14-pl5321hdfd78af_3
  t-coffee           bioconda/osx-64::t-coffee-11.00.8cbe486-h26a2512_0
  unzip              pkgs/main/osx-64::unzip-6.0-h3b3adb9_0
  wget               pkgs/main/osx-64::wget-1.21.3-h6dfd666_0

VEP requires gcc and g++ for installation1. Before running the build script, try setting the CC and CXX environment variables to the paths of your gcc and g++ compilers, respectively:
export CC=gcc
export CXX=g++

You'll need to install these compilers if they don't already exist on your system.

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


text:
Pal2nal translation of large multi-fasta files produces a codon translated file where some sequences half length of the average
I did sequence alignment of a large peptide multi-fasta (n= 4991 sequences). The peptide alignment has sequences with the same length and pal2nal went through just fine... except some of the codon sequences are at half length. If average is X then some sequences are X/2. This is choking IQ-Tree.
I have tried both MUSCLE super5 and MAFFT. The error remain the same (i.e. MUSCLE or MAFFT both lead to some sequences having half of average length) except for different average lengths and average half length in MUSCLE and MAFFT codon sequences. I have pulled out and played with the sequences causing the issue and they seem to be in frame.
Example of peptide sequence not causing an issue: RKVEAFLLFKEMGERGCQPNVHTYTVLIDSFCKERNLDDARKLFDDMFKKGLVPSVVTYNALIDGYCKEGMTEAALEILGMMESKKCNPNARTYNELICGFCKAK
corresponding cds
AGGAAAGTGGAAGCTTTTCTACTTTTTAAAGAAATGGGTGAAAGAGGTTGTCAGCCTAATGTTCATACATACACTGTGCTTATTGATTCCTTCTGTAAGGAAAGGAATCTTGATGATGCCAGGAAATTGTTTGATGACATGTTTAAGAAAGGTTTGGTTCCCAGTGTGGTCACTTATAATGCTTTAATTGATGGGTATTGTAAAGAGGGAATGACTGAAGCTGCATTAGAAATTTTAGGTATGATGGAATCAAAGAAATGCAACCCTAATGCTCGGACCTACAATGAATTGATCTGTGGATTTTGTAAAGCTAAA'
Example of peptide causing issue:
GLCKGGRLNDAWEIFQYLLAKGYQLNVHTYNAMVHGFCKEGLLDEAISLLYKMEENGCVPNSVTFNVVL
corresponding cds
GGTTTGTGCAAAGGTGGTAGATTAAATGATGCGTGGGAGATTTTTCAGTATCTTTTAGCGAAAGGTTATCAACTAAATGTCCATACATATAATGCGATGGTTCATGGTTTTTGCAAAGAAGGTTTGCTTGATGAAGCAATCTCCCTGCTTTATAAAATGGAAGAGAATGGTTGTGTCCCTAATTCTGTAACTTTTAATGTAGTCCTT
Any idea what might be going on?

The example looks fine to me. The peptide causing an issue is 69 residues in length and the cds is 207, which is correct. If the peptide is supposed to be ~140 residues we'd need to see the full length peptide and/or nucleotide. BTW what is the proportion of half lengths against full transcripts?

The peptides are equal in length once they are aligned alongside their orthologs/homologs (i.e. on the other side of sequence alignment). With MAFFT the peptides become something like 31k (including gaps) characters long. Post translation the sequences are something like 93k long and the ones out of whack are about 54k (all of them, the same exact length). Sorry for the confusion, the peptides are not supposed to be the same length per-alignment.

I can see whats happening that the 'issue' sequence is a locus within the full gene for Acacia spp.. The thing is I can't see the original sequences, i.e. before Pal2Nal. The reality I can probably solve this with all the information. At present I can make a guess, but it's no better than a guess.

I am happy to provide the original data, it's all open-source. Do you have a preferred way to receive them?

I don't want to wade through 5K sequences. What I am really asking for is a reproducible example. The example provided doesn't help understand the issue, it's the input into pal2nal I'm looking for. What I want to do is compare the input vs the output.

The first 8 sequences, containing both 'working' vs 'issue' sequences are here. It's a fork of Firefox send, I can use something else if you are uncomfortable with that. By locus if you mean alternate splicing then it is possible that that is what is at play.

Not a problem, doing global alignment on a phylogeny that spans a billion-years was a gamble to begin with. I am grateful you took the time to look at it. I suspect these homo-logs are not real or if they are then the global alignment algorithms we have might not be ready to work on something a phylogeny this deep.

I had a quick blast of the example sequences and this is the Acacia so ancient vascular plant. It forms a locus within the protein in the first part of the question. This protein is comprises a large repeat amino acid motif (check NCBI), the second "issue sequence" described in the question was close to the size of the repeat locus, but start/stop of the locus is not in accordance with the repeat structure. Thus I don't think the repeat locus is part of the issue.
Genome assemblies I noted the examples you gave were all draft genomes and this suggest the problem was with the original genome assemblies. If this was true the way to overcome this would be to reassemble the sequence using a reference genome, bwa-mem2 type thing. Thus, here the error is occurring way before the sequences is even downloaded. Just to note, repeat structures can often be collapsed on a draft genome assembly, but I don't think thats happened here.
Reading frame what I think you are doing is performing an alignment in the proteins and then using that alignment as the template which is then transferred to the sequence data. So now the nucleotide data is nicely aligned including all indels. That appears to be the goal of pal2nal. A common error here is reading frame error, particularly at the beginning of the sequence. I use transalign regularly for this purpose and do encounter that error. I don't know how pal2nal deals with reading frame errors so I can't comment, but thats why I wanted to see the nucleotide input to pal2nal to check this.
Trees the reason this is taking place is because iqtree doesn't permit a tree. That I can help with, I know about tree building. If you post your iqtree command as a separate question , together with the iqtree error, I can probably answer that. What I'll identify is how to build a tree, this likely to exclude sequence, or could be addressed as two trees. Dunno, but if the bottleneck is can't build a tree - thats easily fixable. In this case, is the sequence issue an assembly error, reading frame error ... something else. There's not enough data that I could point to as the underlying source.


text:
Is it possible to do homology inference across species using different kinds of NGS data?
Background: I have a list of species that I want to put through homology inference. The goal of homology inference is to investigate the evolution of a trait on a species tree. I want to use the quantitative aspect of the data like CNV and recreate nodal CNV using cafe or count. And use gene trees of the homology to see if I can discern divergence patterns across gene families of interest.
The issue is some species only have RNA-seq generated transcriptions (peptides and cDNA) while others have genomes. The issue of missing introns aside, how can I circumvent this gap?
Should I:

Just use the cDNA data with genomic data? compare nucleic acid to nucleic acid. I need to eventually make gene trees as well and as such MSA is important.The list of species is huge and I anticipate seg-faults during MSA if not during homology inference itself.

Use emboss or something else to convert the nucleic acid data to peptides.And then compare peptides to peptides? Which makes the homology inference → MSA (hopefully pal2nal and gene tree inference later) a little easier. I have never done translation before so I do not know the pitfalls per se.

My list of species is green plants.
Can anyone share their experience if you have done something like this before?

I think you need to use proteins/amino acids.

Personally, I would use blastx against a local protein database. This converts query cDNA->amino acids and therefore combines points 1 and 2 without need to perform a translation in advance.
Genetic distance and saturation might become an issue here depending whether the search is outside the family and especially Order

MAFFT vs Muscle Super5. I would try "Muscle5" or rather the latest Muscle under the "Super5" algorithm. It outperforms everything.

Thanks for suggesting Muscle5, I will have to try it out this time.

this is the right workflow. deleting my answer.

This, or tblastx to compare translated peptide which avoids the need for a local protein database.


text:
Active sites of the query protein does not align with the top hit using HHpred
I am having difficulties with an alignment on HHpred and was wondering if anyone had any answer to the following:
I am trying to align my protein (PDB 2VU9, a Botulinum Neurotoxin Type a) with one of the top hits generated from Vast.
I know that both my query protein and the other protein both have Glutamic acid as their active site. But, my protein has its active site on the 234th amino acid, whereas the other protein has the active site on its 230th amino acid.
So my question is:
What would be the reasoning for which the top matches on Vast do not align at the active site with my protein, despite sharing the same amino acid at its active site?
Any response helps!

Abstract
"HHpred is a fast server for remote protein homology detection and structure prediction and is the first to implement pairwise comparison of profile hidden Markov models (HMMs). It allows to search a wide choice of databases, such as the PDB, SCOP, Pfam, SMART, COGs and CDD."
Söding et al. Nucleic Acids Res. 2005 Jul 1; 33(Web Server issue): W244–W248.
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC1160169/

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:
Bio.codonalign.codonseq module (cal_dn_ds) for SARS-COV-2
I have calculated the dN,dS ratio for a SARS-COV-2 data set ('NC_045512.2' (Wuhan strain) and 'LC666924.1'), and results in a zero value (for most of the genes).
Is it possible and what is the interpretation, for example negative selection?

Yes that is the expected result, because there are no amino acid mutations in the data sets you are using. Thus the numerator is 0 and anything divided by zero is always zero, even if every codon had a synonymous mutation.
The most precise interpretation is simply there are no amino acid mutations. It doesn't matter how you reword that, all divergence is resulting from synonymous mutations.
The only thing to watch out for is a rounding error, i.e. a single amino acid change is dropping off detection because dN/dS is incorrectly rounded.

A post was raised (now deleted) about what does that say about selection. dN/dS theory states that if it is significantly >1 thats positive selection .. etc. However, a truer representation of the data is there is no amino acid mutation, rather than their is a significant excess of synonymous to non-synonymous mutations. A zero value makes things tricky, and is exemplified by two hypothetical data sets one where every codon had one synoymous mutation (strong negative selection) compared against one that had a single synonymous mutation (little mutation) and both data sets had zero amino acid mutations - both data sets would have a dN/dS of zero but that is not a reflection of negative selection.

To answer the point below a single mutation is not a great basis draw statistical significance and the dN/dS ratio is not a good test if either the numerator or denominator are zero.
Most data is under negative selection in strict accordance with dN/dS.
Neutral selection is extremely rare.
Summary Sometimes the statistical analysis of a data set does not say very much ... and sometimes we have to accept that as an answer.

So if that had a little mutation, Does that mean it's a neutral selection?
And if the sequences (codons and amino) are exactly the same? We can say something about the interpretation?
See examples below-
ORF8
Seq1:
MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVLDFI
Seq2:
MKFLVFLGIITTVAAFHQECSLQSCTQHQPYVVDDPCPIHFYSKWYIRVGARKSAPLIELCVDEAGSKSPIQYIDIGNYTVSCLPFTINCQEPKLGSLVVRCSFYEDFLEYHDVRVVL--I
ORF10
Seq1
MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT
Seq2
MGYINVFAFPFTIYSLLLCRMNSRNYIAQVDVVNFNLT
E
(The same)
Seq1
MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV
Seq2
MYSFVSEETGTLIVNSVLLFLAFVVFLLVTLAILTALRLCAYCCNIVNVSLVKPSFYVYSRVKNLNSSRVPDLLV

I've answered above

@user18071336 If you have a new question do not ask a question as an answer, post a new question, otherwise you will get downvoted and get banned by the system.


text:
Hirschberg's Alignment Algorithm Implementation. Works on the wiki example but not on large sequences of dissimilar size
I implemented Hirschberg's algorithm in python and used the wiki example to verify correct implementation given the scoring parameters and sequences:
Let X = AGTACGCA
    Y = TATGC
Del(X) = -2
Ins(Y) = -2
Sub(x,y) = -1 for mismatch +2 for match

This aligns to:

AGTACGCA
--TATGC-

Which matches the wiki output. I also printed the recursion and it matched the wiki recursive partitioning too. Then I tried it on a harder problem and got a terrible result. It was a mitochondrial Neanderthal DNA sample of about 500 nucleotides and a complete human DNA mitochondrial sample of 16504 nucleotides. The output of trying it on BlastN is here:
https://blast.ncbi.nlm.nih.gov/Blast.cgi#alnHdr_Query_40035
My results look something like this where the top is neanderthal and the bottom is the full human genome. It appears though it just matches each character as soon as it finds on.
---C-CA-----A--------------------G-----T----AT---T--G--A----C-TC--------A--C-C-C-AT--CA-----A--C-----A-----A-C-C-----G-C-CA-T----GT-----ATT--T-C------------------G----TAC----A-T-T-AC---------T----------G------------------C--------CA--------------G-----C--CAC--C-A-T---------GA-AT-AT---------T-----------------------G-T-------ACAG---T-A--C-CAT-------AA--------------------T----------TA-C----TT----------------G-----A--C--T---A-C---C----------T--------------------------------------------------------G----T---A---A---------------T---A-C---ATA-----AA--AACC------------------------T--A---A--T--CC-----A---CA-T-CA---A-----------C---C-C-C--C-CCCC-----C-------------CC-A---T--------GCT-T-A---------CA-A-GCAAGCA-C--------AG------CA------A--TCA--AC---C--------------T--------------T-CA----A---C--T--G--T--C-A---------------------T-A----C----A----T---C-AA-------CTA--C-A-ACT--CC--A---------AA----G----A-C-ACC-C--T-------T-A-C--A---C--------CC-----A---------------C------T----A-----G-------------GA-T---A------TC-A----ACA-AA------C--C-----T------A-C----C----CAC-C-----CT-----------TG--A-------C---A----G-T-A-C---A----T-A---G------CA-CA----T-------A-AA------G---------T-----C--A-------T-----T--T--A---CC-----G-------T-------A--C-ATA------G--CA-C---A-----T-T-----A---TA------G---TC---A--AA----T--------C--C--------C-----T-------T-----C-T-----C--G------CCC---------C----------C-A--T-----------G--------GAT-G-----A-----C--------C--------------C------C-------------C---------C-------T----C----A------G------A---T----A-G-----------GG---------------------G--T--------------C----C----C-T--T--G---A---------------------------------------------------------------------------------------------------------------
GATCACAGGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGATAGCATTGCGAGACGCTGGAGCCGGAGCACCCTATGTCGCAGTATCTGTCTTTGATTCCTGCCTCATCCTATTATTTATCGCACCTACGTTCAATATTACAGGCGAACATACTTACTAAAGTGTGTTAATTAATTAATGCTTGTAGGACATAATAATAACAATTGAATGTCTGCACAGCCACTTTCCACACAGACATCATAACAAAAAATTTCCACCAAACCCCCCCTCCCCCGCTTCTGGCCACAGCACTTAAACACATCTCTGCCAAACCCCAAAAACAAAGAACCCTAACACCAGCCTAACCAGATTTCAAATTTTATCTTTTGGCGGTATGCACTTTTAACAGTCACCCCCCAACTAACACATTATTTTCCCCTCCCACTCCCATACTACTAATCTCATCAATACAACCCCCGCCCATCCTACCCAGCACACACACACCGCTGCTAACCCCATACCCCGAACCAACCAAACCCCAAAGACACCCCCCACAGTTTATGTAGCTTACCTCCTCAAAGCAATACACTGAAAATGTTTAGACGGGCTCACATCACCCCATAAACAAATAGGTTTGGTCCTAGCCTTTCTATTAGCTCTTAGTAAGATTACACATGCAAGCATCCCCGTTCCAGTGAGTTCACCCTCTAAATCACCACGATCAAAAGGAACAAGCATCAAGCACGCAGCAATGCAGCTCAAAACGCTTAGCCTAGCCACACCCCCACGGGAAACAGCAGTGATTAACCTTTAGCAATAAACGAAAGTTTAACTAAGCTATACTAACCCCAGGGTTGGTCAATTTCGTGCCAGCCACCGCGGTCACACGATTAACCCAAGTCAATAGAAGCCGGCGTAAAGAGTGTTTTAGATCACCCCCTCCCCAATAAAGCTAAAACTCACCTGAGTTGTAAAAAACTCCAGTTGACACAAAATAGACTACGAAAGTGGCTTTAACATATCTGAACACACAATAGCTAAGACCCAAACTGGGATTAGATACCCCACTATGCTTAGCCCTAAACCTCAACAGTTAAATCAACAAAACTGCTCGCCAGAACACTACGAGCCACAGCTTAAAACTCAAAGGACCTGGCGGTGCTTCATATCCCTCTAGAGGAGCCTGTTCTGTAATCGATAAACCCCGATCAACCTCACCACCTCTTGCTCAGCCTATATACCGCCATCTTCAGCAAACCCTGATGAAGGCTACAAAGTAAGCGCAAGTACCCACGTAAAGACGTTAGGTCAAGGTGTAGCCCATGAGGTGGCAAGAAATGGGCTACATTTTCTACCCCAGAAAACTACGATAGCCCTTATGAAACTTAAGGGTCGAAGGTGGATTTAGCAGTAAACTAAGAGTAGAGTGCTTAGTTGAACAGGGCCCTGAAGCGCGTACACACCGCCCGTCACCCTCCTCAAGTATACTTCAAAGGACATTTAACTAAAACCCCTACGCATTTATATAGAGGAGACAAGTCGTAACATGGTAAGTGTACTGGAAAGTGCACTTGGACGAACCAGAGTGTAGCTTAACACAAAGCACCCAACTTACACTTAGGAGATTTCAACTTAACTTGACCGCTCTGAGCTAAACCTAGCCCCAAAC

The neanderthal sample matches at the end of the human mitochondrial dna. I suspect the gap score penalty is far too costly to insert enough gaps to reach the end of the alignment but I set blast to a linear gap penalty. How do I overcome this? Here is the implementation code for my program:
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

# NWScore
def NWScore(X, Y):
    # Initialize NM Scoring Matrix. Array needs to be the length of X +1 and Y + 1
    score = [[0 for i in range(len(Y)+1)] for j in range(len(X)+1)]
    for j in range(1,len(Y)+1):
        # Insertion scoring penalty column
        score[0][j] = score[0][j-1] + ins("a")
    for i in range(1,len(X)+1):

        # Delete Penalty Row
        score[i][0] = score[i-1][0] + dele("a")
    for i in range(1,len(X)+1):
        for j in range(1,len(Y)+1):
            # Scoring Matrix
            scoreSub = score[i-1][j - 1] + sub(X[i-1], Y[j-1])
            scoreDel = score[i-1][j] + dele(X[i-1])
            scoreIns = score[i][j - 1] + ins(Y[j-1])
            score[i][j] = max(scoreSub, scoreDel, scoreIns)
    if toPrint:
        print(X,Y)
        print("Score")
        for row in score:
            print(row)
    lastLine = score[-1]
    return lastLine

# Hirschberg
def Hirschberg(X, Y):
    Z = ""
    W = ""
    # If sequence X has no more char in the seq
    if len(X) == 0:
        for i in range(len(Y)):
            # Concatenate a gap(?) to Z and Y[i] to W
            Z = Z + '-'
            W = W + Y[i]
    # If sequence Y has no more char in the seq
    elif len(Y) == 0:
        for i in range(len(X)):
            Z = Z + X[i]
            W = W + '-'
    elif len(X) == 1 or len(Y) == 1:
        #If one seq has a single use full global alignment. For now there is a biopython implementation of global alignment

        a = pairwise2.align.globalms(X,Y,2,-1,-2,-2)
        Z = a[0].seqA
        W = a[0].seqB

    else:
        xlen = int(len(X))
        xmid = round(len(X)/2)
        ylen = int(len(Y))
        # Ok so this is the meat of the algorithm. The last few conditionals dealt with the base cases.
        # Scoreleft takes NWScore for the first half of sequence X, and Sequence Y

        ScoreL = NWScore(X[:xmid], Y)
        # I imagine this is to start from the right side.
        ScoreR = NWScore(X[xmid:][::-1], Y[::-1])[::-1]

        ymid = numpy.argmax([x+y for (x,y) in zip(ScoreL,ScoreR)])

        # Recursion. This is divide and conquer. Creating partitions that somehow dynamically reassemble even
        # though portions of it are getting reversed?!
        Z = Z+ Hirschberg(X[:xmid],Y[:ymid])[0]+Hirschberg(X[xmid:xlen],Y[ymid:ylen])[0]
        W = W+Hirschberg(X[:xmid],Y[:ymid])[1] +  Hirschberg(X[xmid:xlen],Y[ymid:ylen])[1]

    return Z,W
seq1 = SeqIO.read("seq1.fasta","fasta")
seq2 = SeqIO.read("seq2.fasta","fasta")

Z,W = Hirschberg(seq1.seq,seq2.seq)
print(Z)
print(W)
score = 0
for i in range(len(Z)):
    if Z[i] == '-' or W[i] == '-':
        score+= insertionPenalty
    elif Z[i]  == W[i]:
        score+=match
    else:
        score+=mismatch

print("Score: ",score)

The good thing is, this sure ain't a contaminant :-) Is this the untranslated locus of mtDNA?

I believe so. The fasta file states that it is a complete human mitochondrial dna sequence.

Okay thanks, aligning via amino acids will not work then (which is what I woudl have done). The coding looks fine and its a nicely focused question.

Comparisons First of all I would simply say pipe it through Muscle and then compare the results.
I forget the algorithm advances in the 'new' (now quite old) generation of aligners, but it is worth reading up on the latest Clustal omega, possibly MAFFT and how they view Muscle. I do remember is 2018 Clustal Omega upgrade eulogising the historic innovation of Muscle around accuracy and describing how this was implemented. If that works for you, you simply replicate it in your code.
The modern development in aligning algorithms is how they compare in speed versus number of sequences, but accuracy is always measured.
Amino acids -> nucleotides The reason I don't follow the 'accuracy' issues is because its easier to align via amino acids and then switch back to codons. If you are sticking with this gap, penalty system using amino acids is a good approach and then simply back-translating, but if its a UTR then that will not work.
Brute-forcing the gap penalties is certainly an option (which is the advantages of coding it). Proving why it has worked might be a bit tricky, but its doable. You could simply fire parameters onto the data and use regex (import re) or straight character matching == to assess the improvements. Regex is better to assess improvements, whilst == approach is much faster. It depends on the parameter space (I think regex). Moving to an coded 'optical analysis' is an easy modification and could be useful.

The implication from the comment is there's suboptimal code ...

The raw default output to Biopython's pairwise should be used to debug this. Then parameters changed stepwise to understand the anomalous output.

Skimming the code I don't spot anything obvious. What I don't understand is why the code has been duplicated
from Bio import SeqIO
from Bio import Align
from Bio import pairwise2
insertionPenalty = -3
ambigSub = 0
match = 5
deletePenalty = -3
mismatch = -1

#Ambiguous nucleotide meaning
ambNucleotideDict = {"N":["A","C","G","T","U"],
                     "R":["A","G"],
                     "Y":["T","C"],
                     "K":["G","T"],
                     "M":["A","C"],
                     "S":["G","C"],
                     "W":["A","T"],
                     "B":["C","G","T"],
                     "D":["A","G","T"],
                     "H":["A","C","T"],
                     "V":["A","C","G"]}
toPrint = False
#Scoring functions
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

def sub(x, y):
    if x == y:
        return match
    elif x in ambNucleotideDict:
        if y in ambNucleotideDict[x]:
            return ambigSub
    elif y in ambNucleotideDict:
        if x in ambNucleotideDict[y]:
            return ambigSub
    else:
        return mismatch

This is present twice. I think Python would complain about the duplication of functions.
The duplication of variables is fine and probably imports.

A looked through the code a final time and yes there are bugs
def ins(y):
    return insertionPenalty

def dele(x):
    return deletePenalty

That can't be right because whatever you call the functions (which you do), each time the output is -3
No matter what you input the output is always -3. If this was OOP and these were class methods in some sort of method chain that  might be different. These are functions which are operating not different to calling the variable direct. Thus if you substitute
... + dele("a")

with
... + deletePenalty

or even
.... + in("a") # I know this isn't supposed to be here but its to make a point

The answer in all cases is -3.
If you asking for bugs, thats got to be involved.
Out of interest are you using an IDE with breakpoints etc...?

I know there are better ways to do it but I successfully used blastn to align the sequences which operates using scoring and global alignment. So it should be possible. I believe I need to create a system to allow for gaps at the beginning and the end of the smaller sequence.

Okay, the implication is there is a bug in the code, I had a brief look.

My mistake, I removed the duplicated code.

Okay, I looked through the code again and in and dele would a key source of bugs

Large sequences of dissimilar size? Sounds like you want a local alignment algorithm. Hirschberg is a global alignment algorithm, which will often perform unexpectedly on large sequences of dissimilar size.


text:
Can Primer-BLAST be run locally?
Primer-BLAST is a tool for designing target-specific primers, available through the Primer-BLAST web interface.
Is it possible to run the tool locally instead of accessing it via the internet?

Primer-BLAST obviously uses NCBI Blast which we all know about, but it is a pipeline running through Primer3. This is an extremely popular and reliable primer design software, but only does single stranded DNA/RNA. There is a prettier version, Primer3plus here.

To strictly answer the question its not possible because it is a pipeline based on Primer3 which is not written by NCBI, therefore they would not be in a position to distribute it. Nor would Primer3 be able to distribute Blast because that is NCBI. It is simply a web-server.
NCBI and Primer3 would need to come to a joint agreement, which they must have done to present the webserver (at least it would be polite). You can of course replicate the pipeline yourself.
Brief details are here,
https://www.ncbi.nlm.nih.gov/tools/primer-blast/primerinfo.html

It can be done via Docker. This exact pipeline hasn't been distributed via Docker but Blast is available as a Docker image.

New Answer
Just to follow this up, primers4clades is a Docker image aimed at very similar objectives to Primer-Blast. I have not used it, but its worth considering if standalone is essential, e.g. for ease of building a pipeline.
Here's the blurb

Primers4clades provides a fully automated pipeline for the design of PCR primers for cross-species amplification of novel sequences from metagenomic DNA or from uncharacterized organisms belonging to user-specified phylogenetic lineages. It implements an extended CODEHOP strategy based on both DNA and protein multiple alignments of coding genes and evaluates thermodynamic properties of the oligonucleotide pairs, as well as the phylogenetic information content of predicted amplicons, computed from the branch support values of maximum likelihood phylogenies.

I don't know the user, but Docker is a safe running environment that is isolated from the system. There's no harm in checking.
Finally both Blast and Primer3 are available via conda, but not Primer-Blast. In theory, I believe Primer-Blast could be distributed via conda under the rules of conda. In fact, I think, but could be wrong, that anyone could rebuild this pipeline and post it back to conda. @merv will know exactly whether this is 'permitable' and is very likely to correct me if I'm wrong.


text:
What can be the reason of getting negative branches lengths after BEAST analysis?
BEAST2 is currently being used for tree reconstruction prior phylogeographic analysis. The sample size and loci are described below.
I thought that BEAST/BEAST2 does not allow negative lengths of branches on prior distribution level so I was surprised when I have looked at my consensus tree using FigTree.

Question What could be the reason of this and how to wrestle with this phenomenon?
I am analysing about 1500 taxa of HAV (subtype: IA, region: VP1) under the next parameters:

GTR+Г,
Strict Clock,
Coalescent constant population.

The tree was generated by TreeAnnotator through
~/beast/bin/treeannotator -heights median -burnin 20 -limit 0.5 VP1_test_bt_ExpPop.trees VP1_test_bt_ExpPop.tree

Is this an issue of the tree that was generated or an issue of the algorithm that makes the layout for plotting the tree?

Not either/or @pallie. The tree drawing program should be excluded.

What I mean is, does the topology of the tree as saved on disk match the topology of the tree drawn in the graph? @M_

@Pallie It is the problem of the tree. For instance, =0.008838086259721895]:0.6701304093024589)[&posterior=0.005774216356351638]:-0.3146123701838377,(((26[&height=7.0,height_95%_HPD=

@M__ 20 is not a number. This is percentage.

I found the answer on ResearchGate from Santiago Sanchez-Ramirez & let me quote it here:

Try using the option -heights ca when you run TreeAnnotator. This stands for "common ancestor trees" and aims at summarizing clade ages across all posterior trees and not only the values for subset of trees that have that clade. If that subset of trees is low in frequency the average node height might end up being older than the direct ancestor.
Here is the explanation for negative branches:
"MCC trees produced by TreeAnnotator can have a descendent node that is older than its direct ancestor (a negative branch length). This may seem like an error but is actually the correct behaviour. The MCC tree is, by default, generated with average node heights across all trees in the sample which contain that clade. The negative branch lengths result when a clade is at low frequency and tends not to occur in those trees that have the MCC tree's ancestral clade (or vice versa). This means the average heights are for the adjacent nodes are derived from different sets of trees and may not have any direct ancestor-descendent relationship."
In this paper you will find more information about the ca option: Looking for trees in the forest: Summary tree from posterior samples

I am avoiding directly answering the question because in summary, I feel, this given issue is akin to project supervision and is simply too much time.
One issue raised however is something called 'burnin'.
Burnin is foundational to MCMC.
In the comments I stated the burnin is too low. The burnin used ...
~/beast/bin/treeannotator -heights median -burnin 20 ...

The OP believes the burnin represents a percentage, i.e. 20%.
Citing the documentation on TreeAnnotator at https://beast.community/treeannotator

BurnIn (as trees)
In this case you just specify the actual number of
trees you wish to discard rather than the number of states.

In this given case, considering this is 1500 taxa, I would estimate the actual burnin  used in comparison to an initial assessment of completion time (technically called 'convergence') to be 0.4% maximum, i.e. well over 0.5 log smaller than the burnin believed to be used.
It is easy to assess empirically e.g. comparing wc infile and wc outfile, where infile and outfield are MCMC treefiles.
A standard burnin is 10% of the total MCMC. The threshold would be independently assessed and if inadequate a higher burnin would be used, of which 20% is reasonable.
To explain the documentation, in an MCMC only a sample of the trees are written to disk (often 1 in 1000, but it can be anything) and the burnin is 10% of the total (unsampled).
Thus factoring the sampling is important in calculating the precise burnin in any MCMC output analysis. Otherwise the burnin that was intended, isn't the burnin used in practice.
Perhaps you can now see my point?  Not the specific technical matter above, but the issue concerning complexity and time commitment.

To address the matter concerning Beast1 vs Beast2 raised as a comment below, i.e. the proposal the original MCMC was performed in Beast2 and Beast2 is exempt. The issue under discussion is TreeAnnotator.
The TreeAnnotator section of the Beast2 documentation is here https://www.beast2.org/treeannotator/

TREEANNOTATOR ....

BurnIn
This option allows you to select the amount of burn-in, i.e., the
number of samples that will be discarded at the start of the run, so
that you are only analysing the part of the trace that is in
equilibrium.

Again the Beast2 documentation is clear, 'the number of samples', not the percentage or proportion. Furthermore note it states 'samples', not 'states', i.e. the number of trees resulting from the sample frequency that was assigned.
A burnin can be performed prior to TreeAnnotator analysis and that is good practice, albeit not essential. This might be where the report output the OP cites was obtained. Furthermore, if this was done here upstream, then no burnin is required for TreeAnnotator. Moreover if a correct upstream burnin was performed, burnin should either be; avoided for TreeAnnotator, or else the total combined burnin calculated and correctly reported. The final output ain't just a tree, there's a lot of important analysis/ 'meta data'.
Note  I agree with the OP's second comment below this answer using information from the comment below the question: a value of -0.31 can't be ignored, i.e. its a notable value for any branch length and very large in the given context.
Its definitely time to move on
OP upvoted, good luck with the project and all the best.

This is correct for BEAST 1, if I memorise this well. In BEAST 2 the parameter -burnin or -b is percentage. Because it is specified in the -help message & you cannot specify any number above 100 without an exception. Moreover, there are numbers of treated & total trees in stdout of treeannotator. For example, from my another analysis with 10% burning: Processing 18669 trees from file after ignoring first 10% = 2074 trees..

So the question is about possible common reasons of such phenomenon.

Dear @Vovin I have addressed the first comment you raised above as an append to previous answer below the separator.

I have addressed the second point you raised in the section called Note

This is not productive for seeking the solution of the question but I think that I should point this out so others do not get confused by your answer. This is certainly percentage & the fact is based on the next evidence: the help message of the programme where it is named as percentage, the log message where the number is equal to my burnin times the total trees number, & the source code where this variable is called as burninPercentage & used in the appropriate way.

Moreover, for instance, treeannotator -heights median -burnin 101 -limit 1 example.trees example.tree gives beast.app.util.Arguments$ArgumentException: Argument, -burnin has a bad integer value: 101. But it works with 99. The example from exeption message: Example: treeannotator -burnin 10 -heights mean test.trees out.txt. It is obvious that this is not a number because it would be too small. I was myself confused in the past by the blog article that you cited but it works as percentage.

Time to move on


text:
What is & how to solve File error: my.xml.state (Remote I/O error)?
I caught the next exception during my phylogeographical analysis in BEAST 2 with GEO_SPHERE. What could be the reason? & how to evade this in the future?
      ...
      856000000     -3662.2647      5969.5577     -9631.8225 42m16s/Msamples
      857000000     -3522.3540      6105.9849     -9628.3389 42m5s/Msamples
      858000000     -3621.7922      6022.9359     -9644.7282 41m54s/Msamples
      859000000     -3463.1050      6107.9766     -9571.0817 41m43s/Msamples
      860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples
File error: VP1_test.xml.state (Remote I/O error)
java.lang.RuntimeException: An error was encounted. Terminating BEAST
        at beast.app.util.ErrorLogHandler.publish(Unknown Source)
        at java.logging/java.util.logging.Logger.log(Logger.java:979)
        at java.logging/java.util.logging.Logger.doLog(Logger.java:1006)
        at java.logging/java.util.logging.Logger.log(Logger.java:1029)
        at java.logging/java.util.logging.Logger.severe(Logger.java:1776)
        at beast.app.beastapp.BeastMain.<init>(Unknown Source)
        at beast.app.beastapp.BeastMain.main(Unknown Source)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke0(Native Method)
        at java.base/jdk.internal.reflect.NativeMethodAccessorImpl.invoke(NativeMethodAccessorImpl.java:62)
        at java.base/jdk.internal.reflect.DelegatingMethodAccessorImpl.invoke(DelegatingMethodAccessorImpl.java:43)
        at java.base/java.lang.reflect.Method.invoke(Method.java:566)
        at beast.app.beastapp.BeastLauncher.run(Unknown Source)
        at beast.app.beastapp.BeastLauncher.main(Unknown Source)

Unfortunately, googling this error with "BEAST" gave me nothing. If without it, there are a lot of irrelevant problems.

Update
I have just found out that this error happened simultaneously onto all my nodes (different processors, RAM, & common hards & SSDs). All logs of BEAST 2 analyses have it. What is stranger, some of the analyses pressed on after this & someones were broken. To the present I got more interested not in the broken ones but could this exception spoil those which continued their work?

You can restart this at the place you left 'cause is Beast2.
  860000000     -3470.6562      6160.7143     -9631.3706 41m35s/Msamples

However, this is almost 1 billion replications and should be sufficient for convergence. The fact it ends on a round number might infer a algorithm limit is in place. If the Beast2 was restarted from this point and it abruptly ended this would be a very likely reason.
The alternative is an error external to the algorithm such as a RAM error - it has been running a long time.
Whatever the reason I would simply perform ESS statistics and convergence diagnostics on what you have now. There's no real reason to continue if it has converged for 10e6 replications maximum and this is much, much bigger than that.

The terminated processes have not spoiled the work again you can simply input the terminated chains into the restart for Beast2. This is because Beast2 uses checkpoints. Brief information on how to resume from the last checkpoint is highlighted here https://www.beast2.org/2014/04/07/checkpointing-tricks.html. Convergence checks are important.
The situation is only is only catastrophic for Beast1 where you would have needed to start from the beginning.

Unfortunately, big number of parameters demands big number of iterations. I have several thousands of leaves, some sequences are the same but with different dates & locations. Moreover, continuous phylogeography with complex region shapes is really a hard task.


text:
Installing modules in perl for Circos
Im trying to install the missing modules of Circos without success.
The error code is:
cpan[1]> install Config::General
Reading '/home/asg28/.cpan/Metadata'
  Database was generated on Mon, 14 Mar 2022 03:29:03 GMT
Running install for module 'Config::General'
Fetching with HTTP::Tiny:
http://www.cpan.org/authors/id/T/TL/TLINDEN/Config-General-2.63.tar.gz
Fetching with HTTP::Tiny:
http://www.cpan.org/authors/id/T/TL/TLINDEN/CHECKSUMS
Checksum for /home/asg28/.cpan/sources/authors/id/T/TL/TLINDEN/Config-General-2.63.tar.gz ok
Scanning cache /home/asg28/.cpan/build for sizes
............................................................................DONE
'YAML' not installed, will not store persistent state
Configuring T/TL/TLINDEN/Config-General-2.63.tar.gz with Makefile.PL
Checking if your kit is complete...
Looks good
Generating a Unix-style Makefile
Writing Makefile for Config::General
Writing MYMETA.yml and MYMETA.json
  TLINDEN/Config-General-2.63.tar.gz
  /usr/bin/perl Makefile.PL INSTALLDIRS=site -- OK
Running make for T/TL/TLINDEN/Config-General-2.63.tar.gz
Can't exec "make": No such file or directory at /usr/share/perl/5.30/CPAN/Distribution.pm line 2225.
  TLINDEN/Config-General-2.63.tar.gz
  make -- NOT OK
  No such file or directory
Failed during this command:
 TLINDEN/Config-General-2.63.tar.gz           : make NO

I'm running Ubuntu in Windows and already tried in another computer

Can you install the build-essential package on WSL (Windows Subsystem for Linux)? Does sudo apt install build-essential work and, if so, does it fix your issue?

@M__'s answer is perfect, but it merits saying that in the 13 years (as of 2022) since its release many other packages in R or Python can generate plots like Circos. I think Plotly calls it chord plot or something, but most call them Circos plots —like a "hoover" is any vacuum cleaner in BE. Some have less functionality but are way easier to use and install...

The simplest solution is to install conda/Anaconda
Then type any one of the following. Start with the top command and if it is unavailable go to the next commend down - only needs work once:
conda install -c bioconda circos
conda install -c bioconda/label/broken circos
conda install -c bioconda/label/cf201901 circos

This is a preferred solution because the make will already be configured along with any other dependency and in addition you can run it directly on Windows without an intervening simulator as described here: https://docs.anaconda.com/anaconda/install/windows/
You could install brew but this would be a bit messy,
/bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"

because you would be using a perl environment within a simulator.
The alternative is to do what @user324810 and install make explicitly and is described here,
https://stackoverflow.com/questions/23293589/perl-module-install-error-cpan
The error is easy enough but there will be future dependencies which will not be installed and a once an for all solution is a better option and diagnosing and installing each dependency in turn.

The question is about Perl. Is there a conda for perl? As far as I know, conda is a python thing although their site does claim it can do other languages.

Yes. Conda is not singly a Python thing, its a complete environment in this case installing independent Perl and a separate environment to prevent it interacting with system Perl. Using 'brew' is a Perl thing but particularly on Mac its nightmare. Conda is the main reason in my opinion there hasn't been an active uptake of Docker

Just to add, trying is believing. Conda is very cool, because it will load up all dependencies and always works, its better than CPAN, albeit CPAN is likely buried in their somewhere. The only issue is uninstalling because the dependencies remain in place, so you really need to create a new environment for each project (which like doesn't always happen) and delete the whole thing to 'uninstall' everything.

"type any one of the following" - this is not a good recommendation. While Anaconda Cloud automatically lists all labels found on package builds, no one should need the alternative labels (outside of someone deliberately reproducing a previous analysis that used one of those builds). Furthermore, Bioconda expects conda-forge to be prioritized, so a proper recommendation should be conda install -c conda-forge -c bioconda circos. Even better, create a dedicated environment.

Good point @merv, (in fairness I did say start at the top), and its certainly a safer recommendation. If its not entirely mirrored, the entire package should be removed IMO . If you have information in this regard on this package it would be good to share it here. Environments are strongly recommended, 100% agreed, but its the first dive into Anaconda it can a lot to take and the OP may not go with the suggestion (I don't think they did) - its just a suggestion after all.

"If its not entirely mirrored" - not sure I follow what you mean. You mean my comment about needing the conda-forge channel? That's in the docs.

Not quite what I meant but thanks @merv I like the docs. Anaconda is too cool and too useful not to be full part of any longer and these are useful.


text:
Phylogenetic tree rooting in shotgun metagenomics
But I have some weeks fighting with this issue about phylogenetic tree building to use in a phyloseq object in order to calculate beta-diversity metrics that takes into account tree distance branches metrics.
I have one tree for Archea and another tree for bacteria (both unrooted), created from IQTREE from the protein alignments obtained by GTDB-Tk software. I tried to use the software Root-Digger, but only works with nucleotide alignments and not protein alignments. Do you know another software capable (working with protein alignments) to infer more objectively the root rather than define arbitrarily a root through ape R package?

Root-digger is described as "a non-reversible Markov model", which would be assumed to be a directional mutation matrix, thus T->A would have a different probability to A->T using Bayesian approach to branch-length estimation
IQ-TREE is an ultra-rapid method for non-parametric bootstrapping which can be essential for extremely large phylogenies.

Hi @MagiBC this really two questions. Could you kindly focus on your main question and submit the second question as a separate request?

Hi @M__, yes sorry are two questions. Ok, thank you and I take into account your suggestion. The second question is at: https://bioinformatics.stackexchange.com/questions/19018/phylogenetic-tree-building-from-progenomes-database-for-shotgun-metagenomics

One option is to use rootstrap from IQTREE version newer than 2.1.3 (http://www.iqtree.org/doc/Rootstrap). However for at the moment I am not capable to run the first toy tutorial example in the abovementioned question. It's not a stable version. I write in google groups of IQTREE. I keep you informed if I had news.

The answer really is that a root is determined using external evidence, rather than from the internal branches.
An algorithm such as root digger can be used to estimate the root placement based on a branch that for example might minimise the anomalous branch lengths in all other taxa. Excessive comparative branch-lengths could indicate an error in rooting the tree and this can be the basis of the estimate. However, it is an estimate.
R libraries such as ape are more likely to simply mid-point root the tree, which isn't really cool, but its easy to calculate.
The three approaches for rooting are either:

Using a known outgroup to the taxa of interest, which unequivocally predates their divergence.
External information such the evolution of key morphological features of the organism.
Published studies.

Point 1 This can be used singly to find the root - a "root-finding tree" - and the tree of interest presented by constructing a second tree without this root in the alignment. For example the Archaea predate bacteria, one tree determines where the Archaea attach to taxa of interest. If the outgroup is very distinct this could cause artefacts.
Point 2 Finding an outgroup much close to the taxa of interest using external information, for example gram-negative versus gram-positive bacteria, e.g. if all the bacteria in your analysis were one or the other and again do the two tree strategy (point 1).
Point 3 Using large kingdom wide phylogenies from published studies to determine where the most ancestral branch in your taxa.
This is not an exhaustive list about how to use external information to determine the location of the root, but it is very widely recognised.

Directional versus reversible matrices
Directional matrix (e.g. most options in Beast) need to know where the root is, or they will 'guess it': ideally and this information can be feed into the algorithm when it is started. This is important if you are doing molecular dating.
In fact, you can do this even with reversible matrices (e.g. RAxML) but it makes no difference to the answer.

Manual re-rooting
If you have identified the root and the phylogeny you have does not represent that root correctly you can "re-root" the tree using the outgroup identified from point 1 to 3. A program such FigTree are cool for re-rooting. Re-rooting is perfectly cool providing you present the methodology that identifies the root taxa or taxon (singular).

Thanks @M_ for your answer that helps me unclog this issue. I will go through the three approaches to find out the best one in my case. Thanks


text:
Using QCTOOL v2 to process UK Biobank .bgen files - why so slow?
I’m currently using QCTOOL v2 to process imputed .bgen files from UK Biobank, however they seem to be processing very slowly. Is this normal?
My command is pretty basic; I’m filtering out a list of SNPs and samples:
/path_to/qctool \
-g /path_to/ukbXXXX_c22_b0_v3.bgen \
-s /path_to/ukbXXXX_c22_b0_v3.sample \
-og /output_path/ukb_c22_filt.bgen \
-os /output_path/ukb_c22_filt.sample \
-excl-rsids /path_to/snps_rem_c22.txt \
-excl-samples /path_to/samples_to_rem.txt

It is currently processing SNPs at a rate of 1.2/s (e.g. 71169/?,57205.8s,1.2/s).
The computational facility I'm using should not limit the speed of an operation.

Have I made any mistakes in my qctool query?
Or any other ideas on why it's running so slowly?

Alternative suggestions on how to process these files would be appreciated for example, although I would prefer to use QCTOOL perhaps I have to use PLINK.

From the website, https://www.well.ox.ac.uk/~gav/qctool_v2/
"QCTOOL is a command-line utility program for manipulation and quality control of gwas datasets and other genome-wide data"

Hi @Cat could you provide information on computer(s) you are using, e.g. are you using a desktop, laptop or the cloud? What's the RAM of the machine, how full is the hard drive and how big are the files, that sort of thing.

crossposted https://www.biostars.org/p/9519364/

Also whether the files (either input files or output files) are on a local or network drive, that can make a huge difference to the speed.

bgen/qctool2 is just pretty slow in general I've found. For performance, convert to pgen format and use plink2, seems to be orders of magnitude faster as long as you have enough RAM to hold the files in memory,

I contacted the QCTOOL team directly and received this response:

The design of BGEN means that when you subset samples the data has to
be recompressed - this is essentially what makes this slow.  (By
contrast you can subset SNPs very quickly without recompression using
bgenix https://code.enkre.net/bgen.)  It is therefore definitely worth
considering not subsetting samples but using a sample
inclusion/exclusion list instead, if your analysis software supports
that.
If you have to subset and want to use QCTOOL - some things to try are:
I typically use the options -bgen-compression zstd -bgen-bits 8 now,
c.f. https://doi.org/10.1101/308296 this is faster but first check
your downstream software supports zstd compression. Use a map/reduce
type pipeline (i.e. chunk data for re-encoding) - this can be
implemented using bgenix and cat-bgen.

I assume you’re working with the latest version?
Have you tried first stripping the SNPs, then stripping the samples?
Ever since working with the imputed data from the UKB I was never
successful in chopping pieces of bgen files using qctool in a timely
manner. If my memory serves me well, qctool manages to strip SNPs
fairly quickly, but is super slow to remove samples. Try PLINK 1.9 or
2.0.
As for “The computational facility I'm using should not limit the speed of an operation.”,
depending on how disk data moves in/out of slave nodes and how busy
the cluster is, I did saw processes to become I/O starved on large
clusters.
Also, consider that most tools that perform association testing can
utilize SNPs  and samples lists. BOLT, SAIGE, regenie, SNPtest… so
there might not be a need to pre-filter. Put “NA” as phenotypes you’d
like to be “removed” from the testing.

In the end I did not pre-filter SNPs or samples - I set samples to NA within the phenotype file, and used a SNP inclusion list in the second stage of SAIGE with the flag idstoIncludeFile.

When subsetting SNPs (and keeping all samples), I find bgenix to be about 100x faster than qctool.


text:
RNAseq SNP discovery: deciding upon filters and dealing with allele expression bias
I am working with non-model plant RNA samples which we have been deep sequenced and analysed using STAR aligner under default parameters.
Aim We would like to conduct SNP discovery of these samples.
Objective Our ultimate goal with this genotypic data is to search for variants (both SNPs and indels) that fall within genes from a cohort of gene families, which may be responsible for phenotypic variation that we have assessed among our sample plants.
Background The RNA alignments use a reference genome of the species I am dealing with, and following GATK best practices and produced a VCF file for each of my samples.
There are three things that give me pause,

Am I correct in merging the .VCF files after calling SNPs, rather than merging .BAM files and then calling SNPs?
What filters are required for MAF, read-depth, base call quality, missing data, snp proximity to each other and proximity to indels, etc.

Should I produce plots to make these decisions?

Lastly how do I deal with specific allele expression skewing the results, my assumption so far is that other studies seem to have accounted for this by merely using a low bar for the MAF filter (such as >= 0.1), thus ensuring that as many true SNPs as possible are captured.

Any insight into any of these questions, or a point in the direction of where I might find it, would be much appreciated.


text:
How do I set my path for the SRA toolkit under macOS 12.1?
I would like to add the SRA toolkit to my path under macOS 12.1.
After downloading the SRA directory, decompressing, and adding to the path, the directory does not seem to be available e.g. I'm attempting to add it to my path with
export path=(/Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin)

and it does not work. For example if I type
 `ray@dynamic-xx-xxx-xxx-xx ~ % which prefetch

prefetch not found` macOS gives me prefetch not found.  The result of my path is
`ray@dynamic-10-125-206-18 ~ % echo $PATH
/opt/homebrew/bin:/opt/homebrew/sbin:/usr/local/bin:/usr/bin:/bin:/usr/sbin:/sbin`

Is there a way to add it to my path so it works?

Firstly, you are trying to set a path into the Documents directory? Thats like for documents.

Anyway there could be a permission issue. The quickest way is simply set a path within .bashrc or .profile whatever into ~/bin
and then use
ln -si /Users/ray/Documents/School/bioinformatics/sratoolkit.2.11.2-mac64/bin/prefetch ~/bin

... for each compiled program.
$ prefetch

should then work.
This creates a symbolic link in ~/bin so when you type the program in your commandline prompt it will fire, without wiring it via a custom bin. If there is a permission issue a permission error will be quickly flagged.

The best way is to use 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:
conda install -c bioconda sra-tools
conda install -c bioconda/label/cf201901 sra-tools

If one fails the other one will work. Then simply,
conda activate

SRA is now fully configured and ready to go. To get out of the conda virtual environment,
conda deactivate

If you use Anaconda alot it is better to make a specific environment for SRA and not use the core environment, but thats another question.
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 (but forget), otherwise the commandline version is ~/anadonda3.

I see your using a virtual environment via homebrew. Anaconda is separate, its the same idea, but its better.

Fixed it by installing anaconda. It does look like sra-toolkit is not compiled for Apple Silicon yet so Apple x86_64 miniconda is the only architecture that sra-toolkit will work with on MacOS.

Ah right, the M1 chip compatibiility issue. Good to hear that Anaconda has got around the issue.


text:
Taking too long to download metadata file (estimated size 30MB) from Covid-19 Data Portal
I am working on developing a tool to classify covid variants. So I am trying to download SARS-CoV-2 sequence data along with lineage annotation from Covid-19 Data Portal (I could not find a better source that provides both sequence and lineage annotated data.)
Since it is a very large data (given 4,050,433 samples), my idea is to first download the metadata (.tsv file) with accession ids and respective lineages from the portal, filter them to keep only the variants I am interested in, and pass the accession_ids as a text file that I would like to download to CDP File downloader as suggested in the portal.
However, downloading metadata .tsv file from the portal is taking forever (download speed is in bytes sometimes). I have tried downloading directly as a file, using wget and also using curl. I am not sure if there is a way to download it quicker since it is losing connection in between and restarting the download. Just as a side note, the sequence data is downloading pretty fast at 4+Mbps.
I hope my question is clear. Any suggestions/tips on how can I make this happen is very much appreciated. Thanks in advance.

30MB isn't a lot of data at all, given the issue is the EBI servers and not your institutes provision - you would have noticed this before on other projects if it was an institute issue. It therefore suggests someone else is hogging the EBI capacity and there is no provision to bypass this.
The best solution is contact them directly

ecovid19@ebi.ac.uk

... and state the download size and download speed.
There could be a bot sucking up the database (that is what it sounds like), which would be easy to write given FTP is permitted. They will have the ability to block an IP and NCBI have done this before albeit a long time ago.
Finally filling the user experience form would help EBI which is here. However the issue you described suggests escalation. NCBI are very responsive and I would be surprised if EBI did not respond quickly.
Good luck with the project, sounds interesting.
Generic response below if its just traditional over-demand for the service.

I suspect that the beautiful database EBI has provided (and it is amazing) has not anticipated the level of demand for the database size. They are far from being alone.
Modern web architectures, such as Spotify, dedicate enormous effort to matching demand with service provision, notably using Kubernetes. When the server is overloaded they will simply duplicate their core database onto additional machines, which they will be buying in from AWS/GCP. It can even be done dynamically without intervention. NCBI are certainly heading in that direction.
If the EBI wrote this under Kubernetes (it could be a monolith however), the limiting factor is the budget in computer resource allocation (again you can increase Kubernetes provision dynamically). However, if its a monolith this will not be easy and again ultimately someone has pay, i.e. EBI.
Regardless, the fact you can use FTP means EBI are hosting their own service, GCP doesn't permit this and AWS is almost certainly similar. This would mean they have to physically provided additional in-house machines. I visited EBI some time ago and at that time the then director went to some lengths describing how focused EBI were on matching their in-house storage capacity to exponential rise in sequence data. Sudden popularity of a webserver is not trivial issue if it hasn't been anticipated.

Alternatives
Using the NCBI service COVID service is the alternative if there is an urgent need to get the project moving. You may also wish to use GISAID here.
If you want to use the NCBI service, I provided the detailed instructions here: How to programatically download SARS-CoV-2 fasta from NCBI/Genbank via API?
It is however important to understand the possible underlying causes from EBI. I do recall the then EBI director stating their biggest issue was a user simply switching to Genbank, which from their point of view isn't cool given the amount effort going into service provision.

Thank you for the detailed explanation and alternatives. I have found an alternative to EBI from the answer you have suggested. It was quick to download from NCBI.
I will report/inform EBI about the data speed issues and update their response here.

You are welcome. It is a pity for EBI because I think they try hard, however it is also a learning curve for them

I have gotten a response from EBI saying, "They are aware of these downloading issues with few files and they have allocated my ticket to a colleague and help me resolve it asap"

That's good, fingers crossed. They will be sensitive about you switching to Genbank.


text:
Multiple file Python script
I am currently doing an analysis requiring a multiple-input script. I have a directory filled with 900+ samples.
Each sample is comprised of two-files in a fastq format. This is an example of a sample:
MIP_REV_BC_01_01A_S1_R1_001.fastq.gz
MIP_REV_BC_01_01A_S1_R2_001.fastq.gz

The script has a Read1 & Read2 variable, and I need to make it so the R1 file is inputted for Read1 and the R2 file is inputted for Read2 for every sample. To be clear, the script is supposed to be executed on every sample separately.
This was an attempt:
with open(filename) as file:
    lines=file.readlines()
    for line in lines:
        parts[]=line.split("_")
        if parts[7]="R1":
            r1[].append(line)
        else if parts[7]="R2":
             r2[].append(line)

We were bogged down by syntax errors, but either way, this script doesn't seem to do what we wanted it to. However, it captures the idea of what we're trying to achieve.
Question I want to input a directory into the function so it automatically runs the script on the samples in the directory.
Any idea what I should do?

is the idea that you want to process both files at the same time and not store them in memory?

Edited with explanation below. I don't really understand the objective, because I initially thought you were editing the content of fastaq file, but you seem to be editing a file which contains fastaq filenames. Is that correct? I know the syntax but I don't know really what the goal is.
Note I'm using the pathlib which has a glob function, do consider rglob under pathlib.
import re
import pathlib as Path

# pathlib is super, good to habit to use it
mypath = Path('/Users/fewbugs/workingd')
allfilenames = mypath.glob("*.fastaq.gz")
r1, r2 = [], []
x = re.compile('_')
for filename in allfilenames:
    # the 'r' is super important
    with open(filename, 'r') as file:
    # there is no need for a readlines()
        for line in file:
    # split via re, it is better to compile (faster)
        parts=x.split(line)
    # you must use '=='. REMEMBER lists start at 0 not 1
           if parts[6]=="R1":
    # if you want to append to a list, declare the list first
            r1.append(line)
    # there is no such thing as 'else if' REMEMBER lists start at 0 not 1
            elif parts[6]=="R2":
                 r2.append(line)
    # when testing good idea to trap the stuff you're missing
            else:
                 print ("This one escaped %s" % (line))

Please be CAREFUL with my spacing, everything needs to be multiples of 4. I've not checked the code.

Forget this bit I thought you were processing the content of the fastaq file.
outfile1 = Path(mypath.parent, 'outfile1')
outfile2 = Path(mypath.parent, 'outfile2')

with open (outfile1, 'w') as f:
    for item1 in r1:
        f.write("%s\n" % item1)
with open (outfile2, 'w') as f2:
    for item2 in r2:
        f2.write("%s\n" % item2)

A simpler way?
... it looks like you're just processing filenames...
Might I suggest just globing no need for all this stuff.
mypath = Path('/Users/fewbugs/workingd')
r1filenames = mypath.glob("*_R1_*.fastaq.gz")
r2filenames = mypath.glob("*_R2_*.fastaq.gz")

If that is the aim, thats it! Just be careful about the glob versus rglob.
Idea If you want two paired lists this might be better done in a nested comprehension, e.g. a counter could be inserted at the central glob (*).
That may very well be what you're after. i.e.
R1 = ['stuff_R1_001.fastaq.gz', 'stuff_R1_002.fastaq.gz' 'stuff_R1_003.fastaq.gz']
R2 = ['stuff_R2_001.fastaq.gz', 'stuff_R2_002.fastaq.gz' 'stuff_R2_003.fastaq.gz']


text:
How to identify in a RNA seq data in which sample is in which cell using R using qusage
I am using RNA seq data and have been using QuSage in R in hopes to try to identify different samples and figure out which cell it belongs to.
I am trying to visualize and utilize the results but I am having trouble doing that. Does anyone have an idea?

QuSAGE
Quantitative Set Analysis for Gene Expression method
https://bioconductor.org/packages/release/bioc/html/qusage.html

Please be more specific. What have you tried, and where are you running into problems?


text:
Phasing genotyped data using shapeit: ERROR: Wrong ordering in physical positions
I am to phase my genotyped data with Shapeit (from ped or bed files) and I keep on getting this error message:

ERROR: Wrong ordering in physical positions curr_pos=13915
prev=157299737

Any idea on how to fix this?

From the website:
"SHAPEIT is a fast and accurate method for estimation of haplotypes (aka phasing) from genotype or sequencing data. "
https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html

It appears that your SNPs aren't ordered correctly - can you post some of those lines from the .bim file? Also, it seems like you are using shapeit2. Unless you have a good reason, you should be using the most up to date version, shapeit4.


text:
Selection analysis based on each leaves as foreground with FastCodeML
I am using positive selection analysis for each leaves in a phylogenetic tree.
In the routine calculation via codeml, I will be running model zero for all leaves and then model 2 for each leaves as foreground (marking leaves as #1 in phylogenetic tree), then calculate selection using a likelihood ratio test.
Is there a method for using FastCodeML for each leaves as foreground instead of the branch?

First and foremost, the fastcodeml guys didn't replace codeml via refactorising all the functionality of codeml, they refactorised one particular calculation and hence the confusion. fastcodeml isn't a particularly name helpful in truth, because it is not a faster version of all the functionality codeml, its just one of the calculations.
In the authors defence they (notably) improved the runtime of the most intensive calculation in codeml and thats their point - its proof of concept.
Thus, basically no fastcodeml can not be used in this context. fastcodeml runs the branch-site model of codeml only.
What the OP describes isn't how the branch-sites model works, they are describing the branch model.
The branch-site null model sets dN/dS to 1, and the LRT compares dN/dS likelihood between the designated branch(es) between dN/dS fixed to neutral against allowing dN/dS to be any value. Thus it's measuring whether the branch (branches) deviate from neutrality. What the OP describes via, model 0 one rate for the tree versus model 2 meaning two or more rates for the tree (basically two rates) and that is a different analysis and it is called a branch model.

If fastcodeml can be shoehorned to perform sites and branch-models that would be good to know.

There is an issue about 'tip selection' but the key answer concerns the basic difference between codeml and fastcodeml.