text:
How to compare differential expression of same genes between two groups?
I need a statistically significant "decision" that DEGs from group A are different (or not) from group B (after drug treatment).
For that, I have a mean fold change delta, calculated as the mean of 5 different treated samples minus the mean of the same, but untreated samples.
        A.delta    B.delta
gen1    4.9182385  4.5472662
gen2    3.6354761  3.4474844
gen3    3.7592630  3.2838886
gen4    3.9623264  2.8189166
gen5    3.4742171  3.1231050

In my vision, there is not much difference in fold change delta between A and B for the gen1 for example (4.91 is not very far from 4.54). Same for gen2. While for gen4 there might be something interesting.
What kind of plot may be the most useful for that kind of comparison?

I do feel we've been here twice before with this question https://bioinformatics.stackexchange.com/questions/18732/what-is-the-best-way-to-compare-differences-in-deg-expression-between-two-popula and this one https://bioinformatics.stackexchange.com/questions/18660/what-is-the-best-way-to-compare-treated-control-samples both appear very similar questions concerning mouse models. Originally labelled RED and GRN

If the answers you got before not satisfy you, you can always leave there a comment, or edit your question and clearly indicate what about the problem is that you can't use the proposed solutions.

This data is genuinely paired data and you have stated values are a mean. Thus the raw data can be analysed in your previous questions on the same theme described here:
What is the best way to compare treated/control samples?
And here:
What is the best way to compare differences in DEG expression between two populations?
If this is 'ANOVAerable' (its not a word, but if ANOVA is appropriate) than a two-way ANOVA would provide all the answers (however see below regarding the trend in the data). This analysis will perform a 'A' versus 'B' and each 'gen' compared to the other, as well as the interaction between A/B and 'gen'. You would need the raw data, you can't do anything with means (unless you are assuming the distribution). The distribution within each 'gen' is important (if its normal - conforms to the normal distribution, you're in business).
There is nothing further to add and I do feel we need to draw this series of questions to a close.
A regression model could be explored because there is a clear trend in both data sets, but I'm uncertain it will be significant. Within ANOVA however this might trap the outlier, but that is more complicated (you standardise against the trend thus 'gen1' is comparable to 'gen5' and therein 'gen4' can be identified as an outlier). This was touched on in the answer to the previous question (ols and python)


text:
Wrapping R scripts in WDL and run in docker
I am trying to run an R script as a task in WDL in a docker image. Since this is my first time working with these, I am stuck in something very simple (and in the initial stage) and am hoping somebody here can point me in the right direction.
Here is my dockerfile.
FROM rocker/r-base:latest
RUN Rscript ...r-docker/install_packages.R

Here is my WDL script
version 1.0

workflow Rpipe {
   input {
     String name
   }
  call PrintNameR{
        input:
        name = name
}
  output {
     File outfile = PrintNameR.output_R
  }
}

task PrintNameR {
      input {
           String name
}

command<<<
    R --no-save --args name <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

runtime {
            docker: "firsttrydocker"
       }
output {
     File output_R = stdout()
  }
}

I validated, generated an input json file - edited it to include the string I want "Smith", and ran it.
java -jar ../cromwell/womtool-77.jar validate rtry1.wdl
java -jar ..cromwell/womtool-77.jar inputs rtry1.wdl > rtry1.json
java -jar ..cromwell/cromwell-77.jar run rtry1.wdl --inputs rtry1.json

I get an output:
> args <- commandArgs(trailingOnly = TRUE)
> for (i in args){
+ print(i)}
[1] "name"

I get "name", I wanted to get "Smith" because that's what my input is.
What am I doing wrong? Clearly, it seems that I am not passing the WDL task input to the WDL task command -in the Rscript. What is the proper way to do it?

In order to reference a WDL input within a command block you need to use an expression placeholder. Expression placeholders have the form ~{*wdl expression*}, where the inner part can be an valid WDL expression (ie reference to an input declaration)
When an engine sees a String or a command block, it first looks to see if there are any of these expression placeholders within it. If there are, it will evaluate the inner expression, and replace the entire placeholder ~{...} with the resulting value.
You can see how this can be applied below
String first_name = "John"
String last_name = "Smith"
Int age = 45

String full_name = "~{first_name} ~{last_name} # equals "John Smith"
String name_and_age = "John is ~{age}" # equals "John is 45

# Any valid WDL expression can go within an expression placeholder
String complex = "~{first_name + last_name + age}" # equals "JohnSmith45"

A key thing to note is that WDL command templating and expression evaluation is completely separate from the actual execution of the command. An execution engine will evaluate the command template, swap the expression placeholders with their evaluated values, and  THEN, execute the script.
This is relevant for your use case since you are using heredoc to directly define the rscript you are runnning. If we briefly look at a python example, you can see what I mean:
task hello_world {
  input {
   String name = "john"
  }

  command <<<
   python <<CODE
     # this is valid, because once the placeholder evaluates
     # this command becomes name = "john", which is valid python syntax
     # We did not have to pass ~{name} as an arg to the python command
     name = "~{name}"
     print(name)
   CODE
  >>>
}

Now getting to your problem specifically. You are trying to reference the name input value without using any sort of expression placeholder. To the execution engine, the name value looks like any other string. To swap the name in the command for the value passed into the task, you simply need to wrap the variable name in an expression placeholder: ~{name}.
task PrintNameR {
  input {
    String name
  }

  command<<<
    # I would wrap this in quotes to prevent any side effects of
    # shell expansion
    R --no-save --args '~{name}' <<Rscript
    args <- commandArgs(trailingOnly = TRUE)
    for (i in args){
    print(i)}
    Rscript
  >>>

  runtime {
            docker: "firsttrydocker"
  }

  output {
     File output_R = stdout()
  }
}
```

I figured I had to pass ~{sep=" " name} as an argument
command<<<
R --no-save --args ~{sep=" " name} <<Rscript
args <- commandArgs(trailingOnly = TRUE)
for (i in args){
print(i)}
Rscript

source: Stack Exchange

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

sort: votes


text:
Low mapping reads for RNA sequence
I am new to RNA sequencing analysis. I have RNA samples for bacteria for different treatment. I did contaminant filtration, remove adapters and checked fastQC report: attached below for one sample (fwd and rev) .  After mapping with the reference genome using bbsplit but the mapping is only 0.62%.
Any help/suggestions highly appreciated.
Thanks

Well, when you looked at the sequences of some of your reads, what species did they BLAST to?

In addition to swbarnes2 comment, I'd suggest a more systematic approach.
Kraken2 is a tool for metagenomic analysis of reads, though it can be hard to set up.
mash screen is another method for metagenomic analysis, much easier to set up but noisier.
If you are dealing with contamination, treating your data as metagenomic might help.
Also, as always, ensure that your workflow and reference assembly are correct!!!

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


text:
Database to download with information of bacteria?
I am looking for tidy tabular dataset that contains known bacteria and fungi down to the species level. Basically a table that includes genus and species names and other higher-level information. For example, whether it is gram-positive or negative, or anaerobic etc. A database with API would work too.
I found https://bacteria.ensembl.org/, but this is not so easy to use/automate.

You can try Entrez in Biopython. Although, there is hassle with additional information.
from Bio import Entrez

Entrez.email = "your.letterbox@any.com"
Entrez.api_key = "your_key"

TaxID = Entrez.read(Entrez.esearch(db="taxonomy", term="Hepatovirus A"))["IdList"][0]
Entrez.read(Entrez.efetch(db="taxonomy", id=TaxID))

It is a fair point Entrez has an explicit API as described in detail here https://www.ncbi.nlm.nih.gov/books/NBK25498/


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:
A question on Homer normalization
When using annotatePeaks.pl script from Homer software to create histograms, the output is normalized per bp per peak (on top of normalization to 10 million tags).
What does it mean to normalize per bp per peak?

Hypergeometric Optimization of Motif EnRichment (HOMER v4.11, 10-24-2019), is used for de novo motif discovery and NGS analysis and written in Perl and C++.

Crossposted from here, as there is no reply on biostars.


text:
Nextflow: cannot get file pairs from Channel.fromFilePairs
I'm remarkably still struggling with some basic Nextflow concepts. I'm simply trying to collect all pairs of .cram and .crai files from a directory using Channel.fromFilePairs(...), but I get no output when I try to view() the Channel.
Question: Why does my Channel.fromFilePairs(...) command not return file pairs, but I get a list of files when using Channel.fromPath(...)?
As far as I can tell, I'm following the same pattern from the following examples:

Nextflow fromFilePairs
.bam/.bai example (except without the extra code)

Code
nextflow.enable.dsl=2

workflow {

    /*
     * Channel.fromPath returns the expected files and prints them to
     * terminal.
     */
    Channel.fromPath('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram')
        | view()

    /*
     * Channel.fromFilePairs does not print anything to terminal.
     */
    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}')
        | view()
}

Output (note, no file tuples/pairs)
N E X T F L O W  ~  version 21.04.3
Launching `pairs.nf` [pensive_plateau] - revision: f917a5c41b
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram
/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai
/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram

TLDR, use:
nextflow.enable.dsl=2

workflow {

    Channel.fromFilePairs('/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}')
        | view()
}

Results:
N E X T F L O W  ~  version 20.10.0
Launching `xscript.nf` [special_caravaggio] - revision: a4e53657a0
[A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000194-BL-COL-38470BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000018-BL-COL-46250BL1_vcpa1.1.cram.crai]]
[A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000387-BL-COL-38200BL2_vcpa1.1.cram.crai]]
[A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1, [/path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram, /path/to/adsp/cram/snd10000/A-CUHS-CU000563-BL-COL-47451BL1_vcpa1.1.cram.crai]]

When using the fromFilePairs factory method, the trick is to ensure that the pattern you provide can return a common prefix in order for the files to be grouped as expected. An easy way to check to see if this is working as expected is to directly call the readPrefix helper function. For example:
test_cram = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram')
test_crai = file('/path/to/adsp/cram/snd10000/A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram.crai')

pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.{cram,crai}'

println("My CRAM prefix is: ${Channel.readPrefix(test_cram, pattern)}")
println("My CRAI prefix is: ${Channel.readPrefix(test_crai, pattern)}")

You can see that these results do not produce a common prefix:
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [lonely_jang] - revision: 428df7e55e
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1.cram

But if you change your pattern definition to:
pattern = '/path/to/adsp/cram/snd10000/A-CUHS-CU000*.cram{,.crai}'

We can check that it returns a common prefix:
N E X T F L O W  ~  version 20.10.0
Launching `script.nf` [maniac_volta] - revision: a8acf7b03f
My CRAM prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1
My CRAI prefix is: A-CUHS-CU000169-BL-COL-57413BL1_vcpa1.1

Thanks, @Steve. Very helpful. The Channel.readPrefix(...) is a nice test.


text:
On the same strand, for the same gene, can exons be overlapping?
I want to get a set of exon regions for each protein coding gene.
I extracted a set of relevant information (chromosome, start, end, gene ID, gene name, gene type, exon number and exon ID) from a GTF using:
awk '$3 == "exon"' annotation.gtf \
  | gtf2bed \
  | cut -f1-4,10 \
  | awk 'match($0, /gene_name \"[a-zA-Z0-9\_\-.\/]*\"/, gene_name)
    && match($0, /gene_type \"[^"]*\"/, gene_type)
    && match($0, /exon_number [0-9]*/, exon_number)
    && match($0, /exon_id \"[^"]*\"/, exon_id) {
      print $1,$2,$3,$4,gene_name[0],gene_type[0],exon_number[0],exon_id[0];
    }' \
  | tr -d '"' \
  | tr ' ' $'\t' \
  | cut -f1-4,6,8,10,12 \
  | grep 'protein_coding' \
  > annotation.exons.bed

However, during a sanity check, I noticed that for some genes there are some exons with overlapping coordinates, same gene information and exon number, but different exon ID, e.g.:
chr2    9423650 9424052 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001829522.1
chr2    9423653 9423823 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001830696.2
chr2    9423883 9424960 ENSG00000119203.14      CPSF3   protein_coding  1       ENSE00001932298.1

I did not directly work with exons in the past, but I am aware of alternative splicing.
My understanding was that exons can be skipped during alternative splicing, but I assumed that exon number was unique for bases, within the same gene (on the same strand).
In the example above, the first 2 listed exons have the same exon number, overlapping (but not exactly the same) bases, and different exon ID. Is this due to some additional information that I can retrieve from the GTF, i.e., from:
chr2    9423650 9424052 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000475482.5"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-203"; exon_number 1; exon_id "ENSE00001829522.1"; level 2; protein_id "ENSP00000419744.1"; transcript_support_level "4"; hgnc_id "HGNC:2326"; tag "mRNA_end_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000353331.3";
chr2    9423653 9423823 ENSG00000119203.14      .       +       HAVANA  exon    .       gene_id "ENSG00000119203.14"; transcript_id "ENST00000238112.8"; gene_type "protein_coding"; gene_name "CPSF3"; transcript_type "protein_coding"; transcript_name "CPSF3-201"; exon_number 1; exon_id "ENSE00001830696.2"; level 2; protein_id "ENSP00000238112.3"; transcript_support_level "1"; hgnc_id "HGNC:2326"; tag "basic"; tag "MANE_Select"; tag "appris_principal_1"; tag "CCDS"; ccdsid "CCDS1664.1"; havana_gene "OTTHUMG00000090415.5"; havana_transcript "OTTHUMT00000206843.2";

?
Is this because the different exons belong to different transcripts?
However, shouldn't different exons still be non-overlapping one against each other?

I'm not quite sure what you're asking, but pretty much anything can happen in terms of exons and transcripts:

Two different genes encoded on different strands can share the same genomic region as an exon

Two different genes encoded on the same strand can share the same genomic region as an exon (including identical exon boundaries)

A gene can have multiple transcript isoforms that share the same genomic region on the same frame

A gene can have multiple isoforms that share the same genomic region in different frames

Genes can have variable regions which are modified during somatic development (e.g. T-cell and B-cell receptors)

Gene transcripts can be spliced together in different orders in different isoforms

Exons from one transcript isoform can overlap with exons from another isoform (but not necessarily cover identical regions)

Stop codons can occasionally be "read-through", resulting in a longer protein/peptide than the expected translation

As an example of different isoforms within the same gene exhibiting different translation mechanisms, here are graphical representations of annotated transcript isoforms for Gapdh on the mouse GRCm38 assembly, created in JBrowse. I wouldn't put complete trust in any of these (especially not the ones that have no defined intronic regions), but it's a useful demonstration of the types of things that can happen. Compare, Gapdh-207, where an untranslated region starts halfway through an exon, and Gapdh-209, where the translated protein starts halfway through an untranslated exon in Gapdh-207:

If there is a situation where it seems like exons from the same isoform have overlapping regions, it might be worth looking elsewhere in the genome to see if there's a transcribed gene copy somewhere. That, at least, doesn't seem to be something I'd expect to happen through the normal translation process.
Also, don't expect annotations to be complete. Even on something as [presumably] well-studied as the mouse mitochondrial genome, we can identify clear transcript start/stop points from long polyA cDNA sequences that are completely absent in the genome annotations. This is most obvious for the completely unannotated reverse strand (e.g. transcript starting at around 16,200bp), but also happens on the forward strand as well (e.g. coverage steps within mt-Nd5):

Thank you, I believe my confusion was due to the fact that I was not aware of the exon behavior from your two last points. This I think might explain what it is happening in my case, where the exon 1 from the first transcript shares some bases (but not exactly the same interval) with exon 1 from the second transcript. I used to think that transcripts (on the same strand) might have different exons, but the set of exons from which they could select their subset was composed by non-overlapping exons.

May I know which software or tool you are using to show Gadph exons in the images your provided here

I've updated my answer. I used JBrowse (v1 to be specific, but v2 is the most current version).


text:
LD Clumping vs. Conditional Analysis for independent associations
I am working with some GWAS summary statistics. My ultimate goal is to extract variants to use as instrument variables for Mendelian Randomization analysis, which requires me to find SNPs that are independently associated with the exposure variable of interest.
My current workflow involves extracting SNPs from summary statistics with P-value < 5E-08, then performing LD Clumping with a stringent R^2 cutoff (R^2<0.05) within a locus. Is this sufficient to ensure that my variants are independently associated, or should I also plan on running conditional analysis(e.g. GCTA-COJO) in addition to/in lieu of clumping for my variants?

Please clarify your specific problem or provide additional details to highlight exactly what you need. Input/output examples are really helpful. As it's currently written, it's hard to tell exactly what you're asking.

LD Clumping with $r^2$ filtering seems like a good start, but bear in mind that population structure has a significant effect on association statistics.
I recommend carrying out separate PCAs (or some other dimensional reduction) on your variable groups (e.g. risk/non-risk, or high/medium/low) to make sure that their genomic population structure is similar. Alternatively / equivalently, carry out a whole-group reduction, and confirm that the coordinate centres for each variable group are similar.
If the population structures are not similar (which is very likely), then it's important to [also] exclude structure-specific SNPs from the analysis, as these will confound the results of your analysis. Use a clustering method to identify people with similar genetic structure, and treat the different clusters as different trait groupings for identifying those SNPs.


text:
Funcotator reference file error in GATK4
I am using a function called Funcotator (here) in GATK4. The code I am trying to run is:
gatk Funcotator \
--variant chr21.vcf \
--reference hg38.fa.fai \
--ref-version hg38 \
--data-sources-path funcotator_dataSources.v1.2.20180329 \
--output chr21.funcotated.maf \
--output-file-format MAF

However, it gives me an error message like this:
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/Bam/hg38.fa.fai
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.lambda$getFastaExtension$0(ReferenceSequenceFileFactory.java:253)
at java.util.Optional.orElseGet(Optional.java:267)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getFastaExtension(ReferenceSequenceFileFactory.java:253)
at htsjdk.samtools.reference.ReferenceSequenceFileFactory.getDefaultDictionaryForReferenceSequence(ReferenceSequenceFileFactory.java:223)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.checkFastaPath(CachingIndexedFastaSequenceFile.java:180)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:143)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:125)
at org.broadinstitute.hellbender.utils.fasta.CachingIndexedFastaSequenceFile.<init>(CachingIndexedFastaSequenceFile.java:110)
at org.broadinstitute.hellbender.engine.ReferenceFileSource.<init>(ReferenceFileSource.java:35)
at org.broadinstitute.hellbender.engine.ReferenceDataSource.of(ReferenceDataSource.java:27)
at org.broadinstitute.hellbender.engine.GATKTool.initializeReference(GATKTool.java:439)
at org.broadinstitute.hellbender.engine.GATKTool.onStartup(GATKTool.java:722)
at org.broadinstitute.hellbender.engine.VariantWalker.onStartup(VariantWalker.java:45)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.runTool(CommandLineProgram.java:138)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMainPostParseArgs(CommandLineProgram.java:192)
at org.broadinstitute.hellbender.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:211)
at org.broadinstitute.hellbender.Main.runCommandLineProgram(Main.java:160)
at org.broadinstitute.hellbender.Main.mainEntry(Main.java:203)
at org.broadinstitute.hellbender.Main.main(Main.java:289)

Does anyone know how to solve this problem? I tried to use the reference file downloaded from the GATK cloud, but it still doesn't work and has the same error message.
Also, I wonder if the reference files must be the version provided by the GATK team? The one I am using now is not from the GATK team, and as they force me to provide --data-sources-path, I just copied the one they used on their post. Could that cause any problem in the further analysis?
Thanks a lot ~~~
Here is an update:
I first tried hg38.fa as my reference and got this error:
A USER ERROR has occurred: Fasta dict file file:///Users/data/Bam/hg38.dict for reference file:///Users/data/Bam/hg38.fa does not exist. Please see http://gatkforums.broadinstitute.org/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference for help creating it.

So I generated the hg38.fa.fai, then I got the error message shown above.
Then, I tried the reference file downloaded from the GATK cloud, resources-broad-hg38-v0-Homo_sapiens_assembly38.fasta.fai, but the same error appeared again.
Another update:
I tried to use the hg38.dict as the reference and used the same code:
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.dict \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

However, I got the same error message again:
java.lang.IllegalArgumentException: File is not a supported reference file type: /Users/data/hg38.dict

Are all the files .fai, .dict and .fa in the same directory and have the same names?

Yeah, I think so.

You need to pass the fasta file instead of fa.fai or dict file as reference
    gatk Funcotator \
    --variant chr21.vcf \
    --reference hg38.fa \
    --ref-version hg38 \
    --data-sources-path funcotator_dataSources.v1.2.20180329 \
    --output chr21.funcotated.maf \
    --output-file-format MAF

The directory where you have hg38.fa must contain also the hg38.fa.fai and hg38.dict
Hope this helps
Zhan


text:
How to properly annotate sequence variants and errors in a GenBank file format and how to keep track of successive versions of a GenBank file
EDIT
By thinking over the issue, I realized it was mixing two linked but separate questions:

Q1. How to handle the situation when there are divergences between a GenBank file and an actual or theoretical molecule:

Should the difference be indicated in the feature table? If so, with which feature key? And if not, how should it be indicated or what should be done?

Further, when indicating a divergence using a feature key, which should be considered the reference: the file or the molecule?

Q2. How to keep track of changes between versions of a GenBank file when corrections or updates are made?

I have thus re-titled the question accordingly, and by making further researches into the issues I found what seems to me the appropriate solutions (see my answer below).
ORIGINAL QUESTION
Background I am maintaining a database of sequence files describing a collection of bacterial strains and plasmids that we have in our lab.
I am confronted to two sorts of issues:

some sequence files contain errors (see comments below for examples of errors)
other sequence files describe theoretical sequences of molecules that themselves contain variations (mutations, insertion or deletions) that occurred during the molecular cloning process.

Objective I want to annotate the sequences to mark those errors and variations using the GenBank format. However, from the feature keys of the GenBank format described on the INSDC website, it is not clear to me which feature key should be used in each case, essentially this is either 'variation' or 'misc_description'.
Example To better explain the confusion, here below are what seems to me the most relevant feature keys, and my comment of the problem I have with them:
Feature key           variation
Definition            a related strain contains stable mutations from the same
                      gene (e.g., RFLPs, polymorphisms, etc.) which differ
                      from the presented sequence at this location (and
                      possibly others);
...
Comment               used to describe alleles, RFLP's,and other naturally
                      occurring mutations and  polymorphisms;
                  >>> variability arising as a result of genetic
                      manipulation (e.g. site directed mutagenesis) should
                      be described with the misc_difference feature;

According to the line marked >>> above, the feature key misc_difference should be used...
Feature Key           misc_difference
Definition            feature sequence is different from that presented in the
                      entry and cannot be described by any other difference
                      key (old_sequence, variation, or modified_base);
...
Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                  >>> use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

But now, according to the example marked >>> above, it means the variation (here a deletion) is NOT in the considered entry...
Question So how are variations (mutation, insertion, deletion, errors) correctly described in Genbank given that are actually present in an entry ?

This is a good question. One thing I'm puzzled by are 'some sequence files contain errors' if its an ambiguity that is simply dealt with using the degenerate code, e.g. R for A or G ambiguity, N for total ambiguity. If its formatting errors within the file, I'd personally remove them.

@M__ Thank you for the comment. By "errors" I don't mean "ambiguities", I mean sequence files that are not accurate in the sense that the sequence they contain is different from the sequence carried by the molecule they are supposed to describe as checked by sequencing. One could then just correct or discard the files, but for some of them, for the sake of archiving and keeping track of the developments, I would like to be able to annotate those errors.

@M__ Following up your (temporary ?) answer below: * Examples of errors: (i) a feature is properly annotated, but the sequence given is incorrect, (ii) a feature is correctly mentioned, but the sequence span is incorrect (iii) the name of a feature is incorrect (wrong variant for example). * Examples of "variations": (i) during the cloning process a T mutated to a G (variation = 'mutation') (ii) during the cloning process some bases got inserted/deleted (variation = 'insertion' or 'deletion' or 'indel'). Let me know if more details can help.

@M__ That's true, I can see that... Yes, I think too that the most important aspects are mutations (variations), seconded by errors in the sequence data (over annotation errors). Thank you very much again for looking into this.

okay post-updated, official response below.

Your core question is:

(mutations, insertion or deletions) that occurred during the molecular cloning process.

I asked NCBI directly and their response is:

For natural genetic variation use 'Variation' feature;

For 'artificial' genetic variation (lab-produced by genetic modification) then use 'misc_difference'.

In summary in all scenarios its singly the 'misc_differences' that are used (see Conclusion).
So in my scenario where I want to screen out targeted mutagenesis its the 'misc_difference', i.e. the second option of the two you present. This is also true in your second scenario because it is changes possibly unintentional, possibly random mutagenesis / or just unintentional errors during cloning. Finally, this is also true for SOME of the 'errors' you are reporting in comments:

for the first i, ii and iii that appears to be just correcting stuff - administration, just cross-check with those submitting.
for 'variations' the second i and ii in comments, that falls into the 'misc_difference', because a form of mutagenesis, albeit probably unintentional.

Conclusion
The 'variation' feature description isn't clear given the description you reported against NCBI's official Helpdesk response. Its a fair question, its a clear answer, but doesn't explain why the description is confusing.

On 'errors' of mis-identification, e.g. correctly identified e.g. via local Blast or Blat, e.g. sequenced the wrong target, I would use an internal marker, e.g. returned to contributor for comment. Pass it back to the contributor and ask for their input and depending what they say simply mark the record in the 'Comments section' e.g. originally defined as XX. The contributor could have made all sorts of different errors, e.g. annotation - copied and pasted the wrong stuff - experimental anything and you need to give them a chance to respond. Uploading these records on an external database, i.e. Genbank, needs a joint agreement with all parties concerned because it is certainly valid information for someone.

Note, the official blurb here, doesn't appear informative:

Disrupted genes and gene fragments
Sometimes a genome will have adjacent or nearby genes that seem to be
only part of a protein. In many cases these indicate a possible
problem with the sequence and/or annotation. A related issue is the
presence of internal stop codons in the conceptual translation of a
CDS that looks like it should be a real CDS. These problems may be due
to a variety of reasons, including mutations or sequencing or assembly
artifacts. They can be annotated in a number of ways:
Annotate the gene with /pseudo to indicate that there is a problem
with the gene. Note that this qualifier does NOT mean that the gene is
a pseudogene. (see point 2, below, if it is known that the gene IS a
pseudogene) If multiple gene fragments were present initially, then
add a single gene feature which covers all of the potential coding
regions and add the pseudo qualifier. If known, a note qualifier may
be added indicating why this gene is disrupted, for example:

 1        200     gene
                         gene    phoA
                         gene_desc     alkaline phosphatase
                         locus_tag     OBB_0001
                         pseudo
                         note    nonfunctional due to frameshift

Overview
To properly address the questions, each possible case must be listed:

Case 1: divergence between sequences (sequence in the file vs sequence of the molecule); two cases:
a- (wrong) sequence in the file vs (correct) sequence in the molecule
b- ("correct", theoretical) sequence in the file vs (actual, mutated) sequence in the molecule

Case 2: divergence between a feature annotation in the file vs the correct ("true" or reference) biological feature; three cases:
a- wrong location of the feature
b- wrong feature name
c- variant sequence of the feature

Short answer
Case 1a: The sequence in the GenBank file must be corrected and the version number of the file must be incremented (see Q2).
Case 1b: A new GenBank file consistent with the actual molecule must be created and the misc_difference feature key must be used to annotate the variations where relevant.
Cases 2a and 2b: Feature information in the GenBank file must be corrected and the Modification Date of the file must be updated (see Q2).
Case 2c: The misc_difference feature must be used to annotate the variation and the Modification Date of the file must be updated if necessary (see Q2).
Long answer
Q1. How to handle the situation when there are divergences between a GenBank file and an actual or theoretical molecule?
The misc_difference feature key
First, let's analyse the situation with the misc_difference feature key.
According to the INSDC documentation:
Feature Key           misc_difference

Definition            feature sequence is different from that presented
                      in the entry and cannot be described by any other
                      difference key (old_sequence, variation, or modified_base);

Optional qualifiers   /allele="text"
                      /clone="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /replace="text"
                      /standard_name="text"

Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                      use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

From the Definition above, there is still some ambiguity as to what situation is described. Indeed, "feature sequence is different from that presented in the entry" can describe both of the following cases:

the sequence in the file and/or in the molecule is not the one expected for the feature annotated

the sequence in the molecule is not the one represented in the file

Note: In fact, as we will see below, there should be no divergence between the sequence in the file and the sequence of the molecule.
As for the Comment that says that "the misc_difference feature key should be used to describe variability that arises as a result of genetic manipulation (e.g. site directed mutagenesis)", it does not really clarify which of the two cases above is the one intended. Further, the example given in the Comment, which at first sight seems more enlightening, still carries some ambiguity. It says that to annotate a deletion, one should use:
misc_difference 412..433
                /replace=""

All in all, it implies that the deletion is in the molecule (vs in the represented sequence), that it arose as a result of genetic manipulation and is therefore indicated in the file as explained above.
However, it is very counter-intuitive, and not very practical when searching sequences of variants among sequence files, to have a sequence file that is supposed to describe a molecule and in which the sequence is actually different from the sequence encountered in the molecule it is supposed to describe.
And indeed, this view goes against common observed recent practice where the misc_difference feature key is used to indicate how to go from the sequence represented in the sequence file to another version:

Example 1: https://www.ncbi.nlm.nih.gov/nuccore/D90147.1
     misc_difference 308..310
                     /gene="mucA"
                     /note="conflict: deletion"
                     /citation=[1]
                     /replace=""

The indicated deletion is not in the sequence but is a way to go from the sequence to a previously reported CDS sequence for the annotated gene.
Translation of the CDS sequence in the GenBank file:
MKVDIFESSGASRVHSIPFYLQRISAGFPSPAQGYEKQELNLHEYCVRHPSATYFLRVSGSSMEDGRIHDGDVLVVDRSLTASHGSIVVACIHNEFTVKRLLLRPRPCLMPMNKDFPVYYIDPDNESVEIWGVVTHSLIEHPVCLR*

Previously reported CDS (see: https://www.ncbi.nlm.nih.gov/nuccore/M13388;
hyphen/dash in the sequence added by me to mark the difference)
MKVDIFESSGAS-VHSIPFYLQRISAGFPSPAQGYEKQELNLHEYCVRHPSATYFLRVSGSSMEDGRIHDGDVLVVDRSLTASHGSIVVACIHNEFTVKRLLLRPRPCLMPMNKDFPVYYIDPDNESVEIWGVVTHSLIEHPVCLR*

Ref.: https://pubmed.ncbi.nlm.nih.gov/1902211/

Example 2
In each of the GenBank files below, the misc_difference feature key is used to indicate how to go from the (mutated) sequence in the GenBank file to the sequence of the "wild type".

Escherichia coli uidA gene for beta-glucuronidase mutant E413G, clone M2
https://www.ncbi.nlm.nih.gov/nuccore/AM180563.1
 misc_difference 1238
               /gene="uidA"
               /experiment="experimental evidence, no additional details
               recorded"
               /note="site directed mutagenesis"
               /compare=S69414.1
               /allele="wild type"
               /phenotype="Active GUS"
               /replace="a"

Escherichia coli uidA gene for beta-glucuronidase mutant C464R, clone M4
https://www.ncbi.nlm.nih.gov/nuccore/AM180564.1
misc_difference 1390
                /gene="uidA"
                /experiment="experimental evidence, no additional details
                recorded"
                /note="site directed mutagenesis"
                /compare=S69414.1
                /allele="wild type"
                /phenotype="Active GUS"
                /replace="t"

See also:

Escherichia coli uidA gene for beta-glucuronidase mutant Y468C, clone M6
https://www.ncbi.nlm.nih.gov/nuccore/AM180565.1

Escherichia coli uidA gene for beta-glucuronidase mutant W549R, clone M12
https://www.ncbi.nlm.nih.gov/nuccore/AM180566.1

Example 3: https://www.ncbi.nlm.nih.gov/nuccore/KC999035.4
 misc_difference 137414..137415

                 /note="TT compared with single T in original ISCfr1
                 sequence that would result in two partial transposase
                 proteins. Confirmed as correct in raw sequence data."

Here again the misc_difference feature key is used to describe how to pass from the represented sequence to a previous sequence version of the considered gene.

So the misc_difference feature key is appropriate when the sequence in the sequence file does match the sequence in the molecule but does not (any more) match another reference sequence (possibly referred to by another gene-related feature) because variations in the sequence have occurred due to some genetic manipulation or, as in Example 3, due to discrepancy with previous sequencing results.
Thus, the misc_difference feature key is appropriate in Cases 1b and 2c, and the /replace qualifier should be used to indicate how to get from the represented sequence to the other reference sequence for the designated feature (in Example 3 above, the /note qualifier is used instead to describe the difference).
Other cases
There are three other cases left:

Case 1a: the features in the GenBank file are consistent with the molecule, but the sequence in the file is wrong
Cases 2a and 2b: the sequence in the GenBank file is consistent with the molecule, but information about some features is wrong

Those three cases can all be classified as a GenBank file containing errors (in contrast with case 1b where the "error" is in the obtained molecule and with case 2c where the information in the GenBank file is merely incomplete). A GenBank file containing errors is just not appropriate, and the solution here is simply to correct the GenBank file to accurately describe the corresponding molecule.
All that remains is to answer to the second question about how to track changes...
Q2. How to keep track of changes between versions of a GenBank file when corrections or updates are made?
Versioning and time-stamping GenBank files
GenBank files have two fields that can be used to track their history: the Modification Date field and the VERSION field.
The VERSION field uses an accession.version format, where the accession number is fixed and the version number is incremented only when the sequence data is modified.
See:

https://www.ncbi.nlm.nih.gov/genbank/samplerecord/
https://www.ncbi.nlm.nih.gov/genbank/samplerecord/#ModificationDateB
https://www.ncbi.nlm.nih.gov/genbank/samplerecord/#VersionB

Note that only modifications in the sequence result in a change in the version number. Changes in annotations only affect the timestamp (Modification Date) of the file.
See:

https://www.ncbi.nlm.nih.gov/genbank/sequenceids/

With this information we can now describe how to properly keep track of the successive versions of a GenBank file in each of the cases described in Q1.
Proper way to keep track of modifications in GenBank files
Case 1a: (wrong) sequence in the file vs (correct) sequence in the molecule
The sequence in the file must be corrected and the version number of the file must be incremented by one.
Moreover the respective file versions (the former and the newer) must refer to each other in the COMMENT field of their respective header.
See:

https://www.ncbi.nlm.nih.gov/genbank/release/current/ --> 3.4.2  Entry Organization

Case 1b: ("correct", theoretical) sequence in the file vs (actual, mutated) sequence in the molecule
A new GenBank file and new molecule name must be created to accurately describe the variant molecule. The misc_difference feature key along with the /replace qualifier must be used as described above where it is deemed appropriate.
The gene-related features affected should keep the same /gene qualifier, but an /allele qualifier might be used to indicate the presence of a difference and the /note and /phenotype qualifiers can be used to indicate the reference sequence to which the allele should be compared to and any possible phenotypic consequence of the variation.
The misc_difference feature key can be given the same gene-related qualifier information as the gene where the difference occurs (/gene, /allele, /note...)
In the special case of variants unintentionally obtained during genetic construction or genetic manipulation, an altered name for the molecule and GenBank file might be used to reflect the fact that the corresponding sequence is different from the theoretically expected sequence, and the misc_difference feature key should be used to indicate any difference that occurred during the construction compared to the initially targetted sequence.
Cases 2a, 2b or 2c: wrong or incomplete feature data
The feature annotation in the file must be corrected and the Modification Date field must be updated according to the correction date.
Additional references
Sequence revision history

https://www.ncbi.nlm.nih.gov/genbank/sequencerevisionhistory/

Examples of revision history of some GenBank files

https://www.ncbi.nlm.nih.gov/nuccore/U46667?report=girevhist
https://www.ncbi.nlm.nih.gov/nuccore/NM_029716?report=girevhist
https://www.ncbi.nlm.nih.gov/nuccore/NM_033116?report=girevhist

Added for convenience
Here below is reproduced the description and usage of some feature keys and qualifiers that might be relevant to, or linked with, the question of how to describe and annotate sequence variants.
I grouped those definitions there as it might help one to distinguish the different types of variations and divergences and to understand and chose the appropriate term to annotate sequence variants.
Selected feature keys
gene
Feature Key           gene

Definition            region of biological interest identified as a gene
                      and for which a name has been assigned;

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /function="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /operon="text"
                      /product="text"
                      /pseudo
                      /pseudogene="TYPE"
                      /phenotype="text"
                      /standard_name="text"
                      /trans_splicing

Comment               the gene feature describes the interval of DNA that
                      corresponds to a genetic trait or phenotype; the feature is,
                      by definition, not strictly bound to it's positions at the
                      ends;  it is meant to represent a region where the gene is
                      located.

misc_difference
Feature Key           misc_difference

Definition            feature sequence is different from that presented
                      in the entry and cannot be described by any other
                      difference key (old_sequence, variation, or modified_base);

Optional qualifiers   /allele="text"
                      /clone="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /replace="text"
                      /standard_name="text"

Comment               the misc_difference feature key should be used to
                      describe variability that arises as a result of
                      genetic manipulation (e.g. site directed mutagenesis);
                      use /replace="" to annotate deletion, e.g.
                      misc_difference 412..433
                                      /replace=""

modified_base
Feature Key           modified_base

Definition            the indicated nucleotide is a modified nucleotide and
                      should be substituted for by the indicated molecule
                      (given in the mod_base qualifier value)

Mandatory qualifiers  /mod_base=<modified_base>

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /frequency="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)

Comment               value is limited to the restricted vocabulary for
                      modified base abbreviations;

old_sequence
Feature Key           old_sequence

Definition            the presented sequence revises a previous version of the
                      sequence at this location;

Mandatory qualifiers  /citation=[number]
                      Or
                      /compare=[accession-number.sequence-version]

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /replace="text"

Comment               /replace="" is used to annotate deletion, e.g.
                      old_sequence 12..15
                      /replace=""
                      NOTE: This feature key is not valid in entries/records
                      created from 15-Oct-2007.

variation
Feature Key           variation

Definition            a related strain contains stable mutations from the same
                      gene (e.g., RFLPs, polymorphisms, etc.) which differ
                      from the presented sequence at this location (and
                      possibly others);

Optional qualifiers   /allele="text"
                      /db_xref="<database>:<identifier>"
                      /experiment="[CATEGORY:]text"
                      /frequency="text"
                      /gene="text"
                      /gene_synonym="text"
                      /inference="[CATEGORY:]TYPE[(same species)][:EVIDENCE_BASIS]"
                      /locus_tag="text" (single token)
                      /map="text"
                      /note="text"
                      /old_locus_tag="text" (single token)
                      /phenotype="text"
                      /product="text"
                      /replace="text"
                      /standard_name="text"

Comment               used to describe alleles, RFLP's,and other naturally
                      occurring mutations and  polymorphisms; variability arising
                      as a result of genetic manipulation (e.g. site directed
                      mutagenesis) should be described with the misc_difference
                      feature; use /replace="" to annotate deletion, e.g.
                      variation   4..5
                                  /replace=""

Selected qualifiers
/allele
  Qualifier       /allele=
  Definition      name of the allele for the given gene
  Value format    "text"
  Example         /allele="adh1-1"
  Comment         all gene-related features (exon, CDS etc) for a given
                  gene should share the same /allele qualifier value;
                  the /allele qualifier value must, by definition, be
                  different from the /gene qualifier value; when used with
                  the variation feature key, the allele qualifier value
                  should be that of the variant.

/clone
  Qualifier       /clone=
  Definition      clone from which the sequence was obtained
  Value format    "text"
  Example         /clone="lambda-hIL7.3"
  Comment         not more than one clone should be specified for a given source
                  feature;  to indicate that the sequence was obtained from
                  multiple clones, multiple source features should be given.

/compare
  Qualifier       /compare=
  Definition      Reference details of an existing public INSD entry
                  to which a comparison is made
  Value format    [accession-number.sequence-version]
  Example         /compare=AJ634337.1
  Comment         not to be used for new entries from 14-MAY-2021; This qualifier
                  may be used on the following features: misc_difference, unsure,
                  old_sequence and variation.
                  The feature "old_sequence" must have either a
                  /citation or a /compare qualifier. Multiple /compare
                  qualifiers with different contents are allowed within a
                  single feature.
                  This qualifier is not intended for large-scale annotation
                  of variations, such as SNPs.

/function
  Qualifier       /function=
  Definition      function attributed to a sequence
  Value format    "text"
  Example         function="essential for recognition of cofactor"
  Comment         /function is used when the gene name and/or product name do not
                  convey the function attributable to a sequence.

/gene
  Qualifier       /gene=
  Definition      symbol of the gene corresponding to a sequence region
  Value format    "text"
  Example         /gene="ilvE"

/locus_tag
  Qualifier       /locus_tag=
  Definition      a submitter-supplied, systematic, stable identifier for a gene
                  and its associated features, used for tracking purposes
  Value Format    "text"(single token)
                  but not "<1-5 letters><5-9 digit integer>[.<integer>]"
  Example         /locus_tag="ABC_0022"
                  /locus_tag="A1C_00001"
  Comment         /locus_tag can be used with any feature that /gene can be used
                  with; identical /locus_tag values may be used within an
                  entry/record, but only if the identical /locus_tag values are
                  associated with the same gene; in all other circumstances the
                  /locus_tag value must be unique within that entry/record.
                  Multiple /locus_tag values are not allowed within one feature for
                  entries created after 15-OCT-2004.
                  If a /locus_tag needs to be re-assigned the /old_locus_tag
                  qualifier should be used to store the old value. The /locus_tag
                  value should not be in a format which resembles INSD accession
                  numbers, accession.version, or /protein_id identifiers.

/map
  Qualifier       /map=
  Definition      genomic map position of feature
  Value format    "text"
  Example         /map="8q12-q13"

/note
  Qualifier       /note=
  Definition      any comment or additional information
  Value format    "text"
  Example         /note="This qualifier is equivalent to a comment."

/operon
  Qualifier       /operon=
  Definition      name of the group of contiguous genes transcribed into a
                  single transcript to which that feature belongs.
  Value format    "text"
  Example         /operon="lac"

/phenotype
  Qualifier       /phenotype=
  Definition      phenotype conferred by the feature, where phenotype is defined as
                  a physical, biochemical or behavioural characteristic or set of
                  characteristics
  Value format    "text"
  Example         /phenotype="erythromycin resistance"

/product
  Qualifier       /product=
  Definition      name of the product associated with the feature, e.g. the mRNA of
                  an mRNA feature, the polypeptide of a CDS, the mature peptide of
                  a mat_peptide, etc.
  Value format    "text"
  Example         /product="trypsinogen" (when qualifier appears in CDS feature)
                  /product="trypsin" (when qualifier appears in mat_peptide
                  feature)
                  /product="XYZ neural-specific transcript" (when qualifier appears
                  in mRNA feature)

/pseudo
  Qualifier       /pseudo
  Definition      indicates that this feature is a non-functional version of the
                  element named by the feature key
  Value format    none
  Example         /pseudo
  Comment         The qualifier /pseudo should be used to describe non-functional
                  genes that are not formally described as pseudogenes, e.g. CDS
                  has no translation due to other reasons than pseudogenisation
                  events.
                  Other reasons may include sequencing or assembly errors.
                  In order to annotate pseudogenes the qualifier /pseudogene= must
                  be used indicating the TYPE which can be taken from the INSDC
                  controlled vocabulary for pseudogenes.

/pseudogene
  Qualifier       /pseudogene=
  Definition      indicates that this feature is a pseudogene of the element
                  named by the feature key
  Value format    "TYPE"
                  where TYPE is one of the following:
                  processed, unprocessed, unitary, allelic, unknown

  Example         /pseudogene="processed"
                  /pseudogene="unprocessed"
                  /pseudogene="unitary"
                  /pseudogene="allelic"
                  /pseudogene="unknown"

  Comment         TYPE is a term taken from the INSDC controlled vocabulary for
                  pseudogenes (http://www.insdc.org/documents/pseudogene-
                  qualifier-vocabulary):

                  processed: the pseudogene has arisen by reverse transcription
                  of a mRNA into cDNA, followed by reintegration into the genome.
                  Therefore, it has lost any intron/exon structure, and it might
                  have a pseudo-polyA-tail.

                  unprocessed: the pseudogene has arisen from a copy of the
                  parent gene by duplication followed by accumulation of random
                  mutations.
                  The changes, compared to their functional homolog, include
                  insertions, deletions, premature stop codons, frameshifts and a
                  higher proportion of non-synonymous versus synonymous
                  substitutions.

                  unitary: the pseudogene has no parent. It is the original gene,
                  which is functional is some species but disrupted in some way
                  (indels, mutation, recombination) in another species or strain.

                  allelic: a (unitary) pseudogene that is stable in the
                  population but importantly it has a functional alternative
                  allele also in the population. i.e., one strain may have the
                  gene, another strain may have the pseudogene.
                  MHC haplotypes have allelic pseudogenes.

                  unknown: the submitter does not know the method of
                  pseudogenisation.

/standard_name
  Qualifier       /standard_name=
  Definition      accepted standard name for this feature
  Value format    "text"
  Example         /standard_name="dotted"
  Comment         use /standard_name to give full gene name, but use /gene to
                  give gene symbol (in the above example /gene="Dt").


text:
BEAGLE do not see and/or not use OpenCL
I am conducting phylogeny analyses using BEAST 2 as software tool. It is said that it can be accelerated by BEAGLE-library. I have a personal laptop with Windows & a cluster with CentOS.
I installed both BEAST 2 & BEAGLE onto them. This gives me a boost by 3.6 times. It was said too that using OpenCL with BEAGLE gives me more speed-up. Nonetheless, my calculations takes from twelve hours to 3 days so I got really interested in this. Installing it on Windows were very simple. When I ran beast -beagle_info, I saw:
                        BEAST v2.6.6, 2002-2021
             Bayesian Evolutionary Analysis Sampling Trees
                       Designed and developed by
 Remco Bouckaert, Alexei J. Drummond, Andrew Rambaut & Marc A. Suchard

                   Centre for Computational Evolution
                         University of Auckland
                       r.bouckaert@auckland.ac.nz
                        alexei@cs.auckland.ac.nz

                   Institute of Evolutionary Biology
                        University of Edinburgh
                           a.rambaut@ed.ac.uk

                    David Geffen School of Medicine
                 University of California, Los Angeles
                           msuchard@ucla.edu

                      Downloads, Help & Resources:
                           http://beast2.org/

  Source code distributed under the GNU Lesser General Public License:
                   http://github.com/CompEvol/beast2

                           BEAST developers:
   Alex Alekseyenko, Trevor Bedford, Erik Bloomquist, Joseph Heled,
 Sebastian Hoehna, Denise Kuehnert, Philippe Lemey, Wai Lok Sibon Li,
Gerton Lunter, Sidney Markowitz, Vladimir Minin, Michael Defoin Platel,
          Oliver Pybus, Tim Vaughan, Chieh-Hsi Wu, Walter Xie

                               Thanks to:
          Roald Forsberg, Beth Shapiro and Korbinian Strimmer

--- BEAGLE RESOURCES ---

0 : CPU
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_SSE VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_CPU

1 : NVIDIA GeForce 940MX
    Global memory (MB): 2048
    Clock speed (Ghz): 1.19
    Number of cores: 384
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_CUDA

2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )
    Global memory (MB): 3219
    Clock speed (Ghz): 1.00
    Number of compute units: 24
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_GPU FRAMEWORK_OPENCL

3 : Intel(R) Core(TM) i5-7200U CPU @ 2.50GHz (OpenCL 2.1 (Build 10))
    Global memory (MB): 8067
    Clock speed (Ghz): 2.50
    Number of compute units: 4
    Flags: PRECISION_SINGLE PRECISION_DOUBLE COMPUTATION_SYNCH COMPUTATION_ASYNCH EIGEN_REAL EIGEN_COMPLEX SCALING_MANUAL SCALING_AUTO SCALING_ALWAYS SCALERS_RAW SCALERS_LOG VECTOR_NONE THREADING_NONE PROCESSOR_CPU FRAMEWORK_OPENCL

However, it still uses only resourse #0, even through forced beast -beagle_order 3 -threads 1 filename.xml. Although, in comparison to Linux case, it was more successful.
I have downloaded intel_sdk_for_opencl_applications_2020.3.494.tar.gz from here. Then I unzipped it & run install.sh. Installation was successful. After that, I have decided to go to build folder in beagle-lib & run cmake -DCMAKE_INSTALL_PREFIX:PATH=$HOME .. in order to go on to run make install but I get the next message:
-- JAVA_HOME=
-- JNI_INCLUDE_DIRS=/usr/lib/jvm/java/include;/usr/lib/jvm/java/include/linux;/usr/lib/jvm/java/include
-- JNI_LIBRARIES=/usr/lib/jvm/java/lib/libjawt.so;/usr/lib/jvm/java/lib/server/libjvm.so
-- Not using libtools for plugins
-- Could NOT find OpenCL (missing: OpenCL_LIBRARY OpenCL_INCLUDE_DIR)
CUDA_TOOLKIT_ROOT_DIR not found or specified
-- Could NOT find CUDA (missing: CUDA_TOOLKIT_ROOT_DIR CUDA_NVCC_EXECUTABLE CUDA_INCLUDE_DIRS CUDA_CUDART_LIBRARY)
-- Configuring done
-- Generating done
-- Build files have been written to: /export/home/rinkman/beagle-lib/build

It was solved through manual specification of directories cmake -DOpenCL_INCLUDE_DIR=~/intel/system_studio_2020/opencl/SDK/include/ -DOpenCL_LIBRARY=~/intel/system_studio_2020/opencl/SDK/lib64/libOpenCL.so.1.2 -DCMAKE_INSTALL_PREFIX:PATH=$HOME ... The compiling was successful after that. make test was executed without any problems. I detected the file libhmsbeagle-opencl.so.40.0.0 in ~/lib folder. But beast -beafle_info did not show OpenCL. As there are older versions of BEAST & BEAGLE on the root, I tried to change very roughly startup script to run only my local BEAGLE but nothing changed:
#!/bin/sh

if [ -z "$BEAST" ]; then
    ## resolve links - $0 may be a link to application
    PRG="$0"

    # need this for relative symlinks
    while [ -h "$PRG" ] ; do
        ls=`ls -ld "$PRG"`
        link=`expr "$ls" : '.*-> \(.*\)$'`
        if expr "$link" : '/.*' > /dev/null; then
        PRG="$link"
        else
        PRG="`dirname "$PRG"`/$link"
        fi
    done

    # make it fully qualified
    saveddir=`pwd`
    BEAST0=`dirname "$PRG"`/..
    BEAST=`cd "$BEAST0" && pwd`
    cd "$saveddir"
fi

BEAST_LIB="$BEAST/lib"

if [ -z "$JAVA_HOME" ]; then
  JAVA=java
else
  JAVA="$JAVA_HOME"/bin/java
fi

# use BEAGLE_LIB if the BEAGLE library is not in a standard location
if [ -n "$BEAGLE_LIB" ]; then
    if [ -n "$BEAST_EXTRA_LIBS" ]; then
        BEAST_EXTRA_LIBS="$BEAST_EXTRA_LIBS:$BEAGLE_LIB"
    else
        BEAST_EXTRA_LIBS="$BEAGLE_LIB"
    fi
fi

# Explicitly add /usr/local/lib to library search path to ensure
# beast continues to find beagle when installed here.  (This is
# necessary due to known problems with certain JREs.)
if [ -z "$LD_LIBRARY_PATH" ]; then
    export LD_LIBRARY_PATH=/usr/local/lib
else
    export LD_LIBRARY_PATH="$LD_LIBRARY_PATH":/usr/local/lib
fi

BEAST_EXTRA_LIBS="$HOME/lib/"

# use BEAST_EXTRA_LIBS variable to load BEAGLE and other libraries from non-default locations
# this assumes that the library path contains all these libraries (or are set through LD_LIBRARY_PATH)
#if [ -n "$BEAST_EXTRA_LIBS" ]; then
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#else
#  "$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*
#fi
"$JAVA" -Dlauncher.wait.for.exit=true -Xms256m -Xmx1g -Djava.library.path="$BEAST_EXTRA_LIBS" -Duser.language=en -cp "$BEAST_LIB/launcher.jar" beast.app.beastapp.BeastLauncher $*

Please, could anyone help with this?
P.S. I am a noviciate in unix-systems.

@M__ I installed BEAGLE successfully per se & it does work but I would like to run it with OpenCL.

@M__ Do you mean beast-users?

I better explain this style of parameter rich MCMC takes ages and the issues are how speed this up.

Having given this some thought my suggestion is to install miniconda or anaconda and then install the port to OpenCL, then install beast and beagle.
Anaconda is here: https://www.anaconda.com/products/individual
Then install openCL via conda:
 conda create -n beasty python=3.9
 conda activate beasty
 conda install -c conda-forge pyopencl
 conda install -c bioconda beagle
 conda install -c bioconda beast2

My idea was you'd access OpenCL via python described here. This wouldn't work ... its complicated.
Traditionally OpenCL and Beagle didn't work ... on an Intel Mac

2 : Intel(R) HD Graphics 620 (OpenCL 2.1 )

However, if this is your chip set - I wouldn't give up. In any other scenario  I wouldn't pursue this - personally, below.
One alternative is NVIDA e.g. Cloud computing used via CUDA and Beagle. CUDA will easily install and Beagle stacks straight on top.
Any heavy duty MCMC is very amenable to traditional CPU parallelisation - if you think at it - it would be and thats why GPU - CUDA will work. Personally, looking at your resources I would parallelise under you cluster using a MINIMUM of 8-cores per calculation outside CUDA/OpenCL. If that has not already been done, that is the solution in context in my opinion. Unless you're a Google developer few will have access to a NVIDA cluster in any case sufficient CPU parallelisation should be able to achieve the runtime close to that of a good GPU (its gonna be a lot of CPUs though).
Conclusion
Honestly, my personal advise is that phylogenetics sump clusters no matter what the calculation is being used and that has always been the case. The best solution is to max out on available cores. When eukaryote genome guys complain about 1000 hour runtimes (obviously thats the total CPU parallelised), thats a routine calculation in phylogenetics. Buying in NVIDA doesn't necessarily solve this, because whilst MCMC is massively accelerated for maximum likelihood isn't accelerated to that extent by comparison.
Extra information In case the OP is not aware, the enormous advantage of Beast2 is the ability to resume an MCMC at the point it ended (or some checkpoint therein), you don't need to run the whole thing again. You should seek to exploit that, so a calculation started on a laptop can be finished on a cluster. I can't state sufficiently how powerful that is, it stands to reason because the random number stream is declared and every sample thereafter. Keeping all calculations is standard practice in any case, but especially with a Beast2 run and don't every discard it (hard disk is cheap right?). You might find you have to return to it and its a whole lot better than needing to do the whole thing again.
Its time to move on.

I have gotten `Found conflicts! Looking for incompatible packages.
This can take several minutes.  Press CTRL-C to abor
failed
UnsatisfiableError: The following specifications were found to be incompatible with each other:
Output in format: Requested package -> Available versions
Package libstdcxx-ng conflicts for:
python=3.10 -> libffi[version='>=3.4.2,<3.5.0a0'] -> libstdcxx-ng[version='>=7.3.0|>=7.5.0']`

`beast2 -> beagle-lib -> libstdcxx-ng[version='>=10.3.0|>=9.3.0|>=7.5.0|>=7.3.0']The following specifications were found to be incompatible with your system:

feature:/linux-64::__glibc==2.17=0
python=3.10 -> libgcc-ng[version='>=7.5.0'] -> __glibc[version='>=2.17']

Your installed version is: 2.17`
What does this mean? Which packages should I delete? Of what is 2.17 a version? Please could you help?

I get the same conflict with a spanking new environment :(

beast2 can be installed with beagle & opencl only through adding conda-forge & bioconda channels, & creating new environment from scratch as conda create -n env_name opencl beagle beast2. Nonetheless, conda run beast -beagle_info shows only 1 resource: CPU without OpenCL

There is a way to save states in BEAST 1. Should it be left here as a comment or as the answer correction?

That is certainly interesting @Vovin , no its more than a comment. A separate answer would be cool, because I'm not aware of this and it is important or I can ask a separate question.

In this case, I think I could made a Q&A according to this

Thats perfectly, look forward to this @Vovin


text:
Is it possible to reconstruct an MSA with PSSM?
I am thinking of reconstructing an multiple sequence alignment (MSA) of protein sequences from  a position-Specific Scoring Matrices (PSSM).
Is it possible? I suppose co-evolution information is lost in PSSM so maybe it is difficult?

I edited the question to explain what the acronyms represented.

In short I'm afraid it isn't possible because you lose the linkage per Haplotype which is the core information and vital for any population genetics or phylogenetic analysis.
Position-specific scoring matrix e.g. A .60 C .20 G .20 T 0
I tried to produce a position scoring matrix, but it didn't work out, anyway the allele frequency for each e.g. amino acid are given for each alignment position.
You could produce a consensus sequence which, albeit be hypothetical, it might be of some use to understand the generic behaviour of all the sequences in your sample.
Here's a cartoon describing linkage and its importance,

For a multiple sequence alignment of course AB and ab are amino acid residues within the same locus. PSSM could be used to construct a matrix to analyses an established alignment, but that's a fairly complicated Phylogenetics calculation (and a separate question).

Thank you for the explanation from the perspective of linkage!

You are very welcome @JohnnyTam let me know if you have further questions (posted separately of course as per SE rules).


text:
Trinity de novo assembly not completing
I'm trying to assemble some transcriptomes using Trinity, but am having issues getting Trinity to finish. I've been trying to get this to work for weeks, but am hitting a wall and would very much appreciate any and all help.
I have 2 treatment conditions, with 3 replicates per condition, and they were sequenced on a NovaSeq (using 2x150 base reads). I originally tried assembling all the samples together, but even after running for 8+ days it did not finish.
I am now trying to run just 1 sample from each condition, but even after 4 days Trinity will not finish. It seems like it is stuck at Butterfly, and I'm unsure how to proceed. Some more detail and things I've tried:
These are microbial isolates grown in axenic culture, not metagenomic samples. The microbes are protists though (not bacteria), so their transcriptomes are more complex than prokaryotes, but not as complex as higher eukaryotes like humans/mice.
I do not need to have high sensitivity for lowly expressed transcripts, so I tried setting the "min_kmer_cov" higher. First I set it to 2, which did not help. Then I tried 10 (seemed extreme but I was desperate, but even that has not solved it).
-I am running these on a cluster, and admittedly I am not knowledgeable at all about how to optimally request resources, so I am wondering if that is part of the issue. Our cluster uses slurm, and my resource requests right now are
#SBATCH --mem 100000
#SBATCH -n 24
#SBATCH -N 6".

-I tried running just one replicate from one of my conditions. This did finish (after ~42 hours), so it seems like the issue is not an inherent problem with the transcriptomes/the data itself.
Thank you very much!

By efficiently constructing and analyzing sets of de Bruijn graphs, Trinity fully reconstructs a large fraction of transcripts, including alternatively spliced isoforms and transcripts from recently duplicated genes. Compared with other de novo transcriptome assemblers, Trinity recovers more full-length transcripts across a broad range of expression levels, with a sensitivity similar to methods that rely on genome alignments. This approach provides a unified solution for transcriptome reconstruction in any sample, especially in the absence of a reference genome.

Trinity
Trinity is a de novo assembler that fully reconstructs a large fraction of transcripts, including alternatively spliced isoforms and transcripts from recently duplicated genes using de Bruijn graphs. Trinity recover  full-length transcripts across a broad range of expression levels, with a sensitivity similar to methods that rely on genome alignments.
Grabherr, MG et al. (2011) Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nature Biotechnology. 2011 May 15;29(7):644-52.

i worry a little bit about the co-assembly of multiple replicates. i don't know trinity well and as of 10 years ago it was the recommended strategy, but especially in the presence of splicing or RNA processing funny business (and protists are famous for funny business) you could assemble transcripts that never actually exist. metagenomic assemblers are optimized to avoid this error mode, IDK if trinity is. an alternate approach would be to assemble independently and merge.

I'm not familiar with the use of Trinity for microbial assembly. There may be some assumptions for animal assembly that are broken in protists, although I can't think of any reason why that would slow the butterfly step down.
It's possible that your dataset is simply too large. One thing you could try is reducing the dataset size (after pooling all samples together). Trinity can do a kmer-based normalization, more details here:
https://github.com/trinityrnaseq/trinityrnaseq/wiki/Trinity-Insilico-Normalization

Large RNA-Seq data sets, such as those exceeding 300M pairs, are best suited for in silico normalization prior to running Trinity.
Note, all recent versions of Trinity will perform in silico normalization by default. You can turn it off with Trinity --no_normalize_reads.
If for any reason, you'd like to normalize data separately, you can do so as follows:
 %  $TRINITY_HOME/util/insilico_read_normalization.pl

If that normalization is still failing to reduce the dataset down to a reasonable level, it suggests that the dataset is probably more complex than expected. This could be caused by contamination with other microbes, and may be improved by first filtering out contaminant reads.
In any case, if you're having difficulty with getting Trinity to work, the best approach is to contact the developers by logging an issue on the Trinity github:
https://github.com/trinityrnaseq/trinityrnaseq/issues
You could try using an alternative short-read genomic assembler like SPAdes. I see that SPAdes also has an RNA mode and a metaviral mode, so it might already be a bit more closely tuned to smaller transcriptomes:
https://cab.spbu.ru/software/spades/
https://github.com/ablab/spades/blob/spades_3.15.4/README.md#metapv


text:
Help me to calculate Heaps Alpha value from the roary pangenome pipeline result?
I need to know whether my pan-genome is open or closed. For that, I need to calculate the Heaps Alpha value from my roary pipeline results. Please help me to do the same.
Many thanks in advance.

Directly quoting from the link provide:
Roary By Andrew Page based on version 3.11.2 (22-Jan-2018)
Roary is a high speed stand alone pan genome pipeline, which takes annotated assemblies in GFF3 format (produced by Prokka (Seemann, 2014)) and calculates the pan genome. Using a standard desktop PC, it can analyse datasets with thousands of samples, something which is computationally infeasible with existing methods, without compromising the quality of the results.

I edited the question to give details of Roary. I personally find the context difficult to follow and further details would be helpful, such as the aim of the project as well as the genomes being used.

If this is a homework question (or otherwise canned question with a known answer), could you please tag as 'homework'? Also, please provide additional context around what you're trying to do (e.g. sample type, project goal).


text:
Mapping to assemblies in subdirectories
I am performing a bash loop one liner:.
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

The file hierarchy is as follows:
\Drain_4
    \Assembly
        \A4_assembly
           assembly.fasta
           A4.fastq
        \D3_assembly
        \I8_assembly
        \U7_assembly
\Drain_5
\Drain_7

Where \ denotes a directory. This is the error i am getting:
[ERROR] failed to open file '.fasta': No such file or directory A4.fastq

Both input files are in the same directory.
The fastq is my current directory. I'm guessing i will need to provide the complete pathname for both my input files?
I should mention that in all folders the assembly.fasta file is named the same but the fastq file is different (has a different prefix depending on which folder it is in.

your file hierarchy doesn't indicate any file called NC_D2/NC_D2_trim_filt.fastq; where is this file from?

The most likely issue is that the assembly.fasta file is not found in the current directory. Taking @m's comments into consideration, a minimal rewriting of your script to fix this would be as follows:
for k in */Assembly/*/*.fastq; do minimap2 -ax map-ont $(dirname ${k})/assembly.fasta $k > $(dirname ${k})/mapping.sam; echo $k; done

But I'd recommend at least adding additional quote and brace protection, just in case file names contained spaces:
for k in */Assembly/*/*.fastq;
  fDir = $(dirname ${k});
  do minimap2 -ax map-ont "${fDir}/assembly.fasta" "${k}" > "${fDir}/mapping.sam";
  echo ${k};
done

edit: fixed up variable capitalisation

The solution is quite simple you issue the one liner either from within A4_assembly directory or else NC_D2 (see below). It will then find assembly.fasta and the error will disappear. That's because the current working directory is assigned to the directory in which the one liner is issued. Directory A4_assembly is where assembly.fasta is situated, thus is will work.
If will fail to find the fastq file thus you will need to rewrite the loop, simple,
for k in *.fastq; do minimap2 -ax map-ont assembly.fasta $k > mapping.sam; echo $k; done

... will work, executed from within Assembly_4
The one liner suggests you are executing the command from the directory above the Drain_4 directory.
If this is how you intended the one liner to operate I don't think its the solution to the one liner, because there is a risk you will encounter logical errors (not bugs) ..
*/Assembly/*/*.fastq

If you did execute this commend from above Drain_4 it will looking all over the place here in every Drain and every assembly directory therein. If you executed this anywhere else it will fail to find *.fastq because the start of the wildcards is the current working directory in a bash for relative paths.
However, when you dump the file it goes into one one place:
 > mapping.sam

So there is a risk are performing loads and loads of minimaps but they are simply overwriting each of other, i.e. they are replacing each other.
The other issue is
NC_D2/

This directory is not represented in your directory file.
Fundamentally the problem is that assembly.fasta is always in one location which appears to be in the directory above the directory structure you described, viz. */Assembly/*/*.fastq ... according to the one liner, but in the directory structure its in A4_assembly it is also likely to be in \D3_assembly amongst loads of others and other Drains (master directories).
Thus even if it 'finds' the reference genome (if that is what it is?), it will see that for everything ... and I mean everything.
There are dynamic solutions to what I think the code using bash is looking to achieve because you can issue regex from within the do ... done, or simply do it as a script.

I suggest to simply execute the first script I wrote in every Assembly directory then everything will simply work and there's no risk of logical errors.


text:
Workflow for identification of new splicing variants from RNA-seq data?
I'd like to perform a search for possibly unidentified splicing variants of a specific protein in A. Thaliana. I do not have my own RNA-seq data. It has been suggested to me to use this workflow:

Download RNA-seq data from SRA with fastq-dump.
Construct an index of reference sequences with bowtie2-build.
Map the reads from RNA-seq data to the reference sequences with bowtie2.
Remove non-specific mappings with samtools view -q 10. Sort the mappings with samtools sort.
Browse the aligned reads with tablet.

I have been given example data of another organism (not A. Thaliana) and have successfully completed these steps. However, I don't know how to apply that in my research. Specifically, I'd like to know:

Is there other/better way to do do this?
Can I narrow the search/alignment to a specifc gene?
Where can I download A. Thaliana's reference genome in FASTA format?
Which SRA runs should I choose for my analysis?
How to interpret the results to get the information on the different splicing variants?

Linking a similar question I asked before: How to find an unidentified splicing variant of a protein?

The standard format is to ask one question per Bioinformatics SE question. You might think about focusing one or two questions within the 5 questions you raised.

I'm a begginer in this topic and I'm not exactly sure what I'm looking for. I wanted to show I did some legwork and ask specific questions about the topic instead of one general question "how to do X?". I need to answer all of them in order to progress with my project. What should I have done instead? Create separate posts for each question?

Is there other/better way to do do this?

There is almost always a better way for pretty much anything but in this case I would stick with what your colleagues are using/suggesting. Once you have a worked-out case, you can always improve if needed. Other than that what is suggested to you seems reasonable to me (I have done sth similar for a gene family of some 100 genes using a workflow similar to what was suggested to you)

Can I narrow the search/alignment to a specifc gene?

Technically yes, instead of the whole genome, you can use a subset of the genome. This will present some problems, for example since the aligner will not be able to find native "targets" for most of the reads, some reads will map to your "new reference" (genome subset) in a "False-positive" manner. You will need to play with the settings. I would advise against this.

Where can I download A. Thaliana's reference genome in FASTA format?

Genbank, EMBL, ... There might be a specific repo for your organism like there is Flybase for Drosophila species and Xenbase for Xenopus species.

Which SRA runs should I choose for my analysis?

RNA-seq experiments conducted on your species of interest and ideally on a tissue where your genes of interest are highly expressed. If there is none, you could use a similar organism but this would come with some problems. The most ideal would be RNA-seq with long-read sequencing (like Pacbio).

How to interpret the results to get the information on the different
splicing variants?
When you visualize the aligned reads and compare these with existing gene annotation model, it will immediately be obvious if there are novel isoforms, could be due to "new" or not-yet-known exons/splice signals.


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 can I draw a diagram of hydrogen bonds only on the basis of C-alpha backbone?
We are working on an algorithm to assign proteins only on the basis of alpha-carbons.
My teacher told me to draw a diagram to represent hydrogen bonds only on the basis of alpha-carbons.
I showed him the following picture:

He said this picture doesn't count, as this diagram shows all the atoms.
I couldn't figure out how to do this without showing oxygen, non-alpha-carbon, and nitrogen atoms.
How can I draw a diagram of hydrogen bonds only on the basis of C-alpha backbone?

There are three options:

it's a badly phrased Ramachandran plot question
it's a question aimed at the describing the properties of different SS
it's a trick question

φ & ψ can be used to classify SS (cf. Ramachandran plot), these are however proper dihedrals, so require 4 atoms, which you lack.
The sequential Cαs in a strand are increasingly distant by 3–4Å, while in a helix there's a periodicity every 3–4 residues. The images you were given are likely to remind you that between two strands in a sheet the Cαs are ~5Å away from the other strand, but the order is different —the indices will increase together in parallel strands, but not in anti-parallel strands.
A C-alpha trace has the issue that the residues are often not known (UNK is technically represented as homoalanine, but is an unknown) and the resolution is > 4Å. Any SS allocation will be highly inaccurate making it a trick question.


text:
Constructing Two Triangle Heatmaps in One Square
I'm looking to make a two, two range colour scheme heatmap where each range corresponds to a separate triangle in the heatmap square.
I am looking for a relevant guide for this heatmap. Any recommended package/library or guides? Thanks!

Are you using R or Python, i.e. you need a Dataframe and its what language you are okay with? Secondly, are you fine doing a single colour heat map? For example even I know with ggplot2 this is easy if you're using uniform colour scheme?

I'm okay with both, but I prefer more on Python. It's okay if you want to show me the single colour heat map with ggplot2 (I'm okay with R).

If you're not married to this idea and just wish to show correlations and significance at the same time, R's corrplot can do this. It uses *, ** etc. as content to show p-value thresholds, and a color scale to show correlation value.

This solution doesn't reproduce the plot, but it would do if the numerical ranges between the triangles do not overlap.
import matplotlib.pyplot as plt
import seaborn as sns

xx = plt.get_cmap('rocket').copy()
xx.set_under('none')
yy = plt.get_cmap('mako').copy()
yy.set_over('none')
axis = sns.heatmap(dataframe, vmin=-1, cmap=xx, cbar_kws={'pad': -0.04})
sns.heatmap(dataframe, vmin=0, vmax=1, cmap=yy, ax=axis)
plt.show()

My colour schemes are not great but they are distinct, you can search for more here:
https://seaborn.pydata.org/tutorial/color_palettes.html
The idea is this would address a situation where you compare the r-value OR a p-value against another set of values where are > 1. That's what the vmin and vmax are doing.
This is a ONE np array or dataframe solution. In other words you will need to combine your data sets. The solution you want is a two dataframes in one heatmap square. You can put two triangles side-by-side, i.e. one dataframe per heatmap, but I don't know who to place two separate dataframes in one heatmap.