text:
Plink 1.9 --merge-list changing order of fam file
I am doing some whole-genome sequence analysis and to speed things up split the data into chromosomes and also into regions within those chromosomes. (So I have my dataset split into ~300ish files) I am doing some IBD analysis with plink files and to do this IBD analysis I am combining the files into one big plink file.
I ran into issues today where when using plink 1.9 and the --bfile and --merge-list commands, the output fam file is different from the input fam file.  My dataset has related individuals so the fam file is very important.  Also its is important to note that when using binary files the order in the fam file needs to be exactly the same as the order the individuals appear in the other files (I learned this the hard way haha).
Does anyone know why plink 1.9 changes the order in the fam file when merging plink files, and how to fix this?  It works with plink 1.07 but with 1.9 I get a different fam file.
Performing "diff" on input and output fam files shows they are different:
plink=$mydir/bin/plink_1_9/plink

$plink --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

Using the same code but using plink-1.07 creates matching fam files as expected"
plink1=bin/plink-1.07-x86_64/plink
$plink1 --noweb --bfile $indir/chr22/285/$prefix.chr.22.region.285.qc --merge-list $indir/chr22/merge_list.txt --make-bed --out $indir/chr22/$prefix.chr.22.qc

The merge_list.txt file is in this format:
$indir/chr22/286/$prefix.chr.22.region.286.qc.bed $indir/chr22/286/$prefix.chr.22.region.286.qc.bim $indir/chr22/286/$prefix.chr.22.region.286.qc.fam
$indir/chr22/287/$prefix.chr.22.region.287.qc.bed $indir/chr22/287/$prefix.chr.22.region.287.qc.bim $indir/chr22/287/$prefix.chr.22.region.287.qc.fam
$indir/chr22/288/$prefix.chr.22.region.288.qc.bed $indir/chr22/288/$prefix.chr.22.region.288.qc.bim $indir/chr22/288/$prefix.chr.22.region.288.qc.fam

Can someone explain to me why this could be happening and if there is a way to fix it? Thanks so much!!
(Note, I know I can manually replace the fam file with the correct one with cp or mv, but I would like to know how to handle this with plink 1.9, if possible, as I am afraid it could be indicative of an underlying error in the merge process :) )

Thanks Christopher Chang for the great answer via the plink2 google group! (See here)
His answer was as follows:
"Yes, plink 1.9 did change the default merged sample order.  However, you can request plink 1.07’s behavior by adding “--indiv-sort 0” during the merge.  (You can also use --indiv-sort + --make-bed at any time to change the sample order to whatever you want.)
As for why the default was changed, take a look at the new default order: it will probably make a lot of sense.  It’s “natural-sorted” which tends to make it easy for humans to manually search through regardless of what mix of letters and numbers are used in the IDs."
See also the Plink 1.9 documentation for --indiv-sort which makes the point that this behavior is default only for merges--in other commands Plink 1.9 uses the same ordering as in Plink 1.07, which is the order the samples were in.


text:
Dealing with indels when converting a position to chromStart and chromEnd for a bed file
Starting with a chromosome # and position, I am trying to get chromStart and chromEnd values for a .bed file, but I am not sure how to calculate chromEnd when I have a variant that is an insertion or deletion.
My file is not a VCF file but it gives the chromosome number, chromosomal position (pos), reference (ref) allele, and alternate (alt) allele for each genetic variant, like you would find in a VCF file.  The tricky part is that not all of my variants are single nucleotide base changes--I also have insertions and deletions that can be multiple bases long.
For SNPs I am setting chromStart = pos - 1, and chromEnd = start + 1 = pos.  (Note, .bed format requires 0-based indexing and the end position is exclusive. See https://genome.ucsc.edu/FAQ/FAQformat.html)
I am unsure of how to handle indels though (insertions and deletions).  Currently I am adding the length of the reference allele to the starting position to get the end position (chromEnd = start + len(ref)).  I am not sure if this is right though.
I will be using this .bed file with a .bed file from Gencode and bedtools. Does anyone know the standard way the research community (or at least Gencode) treats indels when determining chromosomal position for .bed files? Thanks!

Usually people don't convert such things to BED files. bedtools understands VCF format, so why not convert to that, since it seems most of the information is already there?

Could you explain why you are trying to make this into a bed file? Is it for displaying it in a genome browser? I don't really see how a bed file could even show an indel, unless it were by using the coordinates of the reference genome.

Hi, @terdon I am trying to make this into a BED file so that I can find which variants overlap exonic areas in the genome. I have a BED file from Gencode listing protein coding regions, and I want to know which of my variants overlap with the exome.

Hi @Devon Ryan, if the best thing to do would be to convert it to VCF format, will the bedtools intersect command be able to compare a VCF file with a BED file, or would I need to convert the Genecode BED to VCF?

bedtools is able to compare BED and VCF files, you don't need to convert anything at that point.

The best way to determine which variants affect exomes is to run through an annotation tool like the VEP (http://www.ensembl.org/info/docs/tools/vep/index.html) for which you need a VCF.

You can convert to VCF and use BEDOPS vcf2bed --insertions, vcf2bed --deletions and vcf2bed --snvs to segregate the VCF into three subsets.
For convenience, BEDOPS conversion tools deal with index shifts, including vcf2bed. The resulting BED file is 0-indexed, half-open.

Thanks @Alex Reynolds.  Is there any way to do this but without segregating the VCF into subsets?

Coming back year later I can say genomic coordinate systems can be tricker than I ever imagined years ago when I wrote this question, especially when converting between different formats lol. Here is a good article explaining 1 based vs 0 based positions. https://www.biostars.org/p/84686/ . VCF format is generally 1 based and BED format is generally 0 based. But also some programs use 1-ish or 0-ish indexing so it gets confusing (looking at you VarSeq ;) )

source: Stack Exchange

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

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

sort: votes


text:
Two variants associated with same chromosomal position in genotype data for one individual
I'm trying to play around with my family's raw genetic data from Ancestry.com so that I can do some genetic analyses, and I noticed in the AncestryDNA raw data files that some regions of the genome have more than one variant associated with them and I would like to know why this is. The end goal is get this into a regularly formatted VCF, but I feel that these double variants will probably mess things up. So for example one might see something similar to this in the tab-delimited file:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

(Header line: rsid    chromosome    position    allele1    allele2)
I understand that when rsid numbers were created, sometimes multiple rsid numbers were accidentally assigned to the same region of the genome, but  I do not understand how two separate variants can be assigned to one single region in the genome. I could be wrong but as I understand it, I should only see one line here (either a "T    T" or an "I    I") because genotypes are reported in comparison to the human reference genome and so each position should uniquely refer to a single nucleotide in the said human reference genome.
But in this example case one line says there is a "T    T" genotype at this position, and the other line says there is an insertion(see "How do I read my raw DNA data?" here).  One possible explanation I thought of is if this would just be saying that a T was inserted, but if this is the case, then as I understand it the reference allele present in the individual should still reported, just in conjunction with the allele at the index. (So assuming an extra "A" was inserted after the T, which could be the reference or alternative allele, something like this should be recorded: rs6010164  22  50515270    TA  TA).
Can anyone help me understand why there are two separate variants associated with chromosomal position? Thanks!

Are you sure that if an A was inserted, you would get GA GA? I have seen G GA for an A insertion. Looking back at your header, it does not seem to respect the VCF format. The order of the first three columns is scrambled. It should be chromosome, position, id. Then, there should be the reference sequence at that position, followed by the alternative sequence. Then, in the sample column, you should get the genotype in function of the reference and alternative alleles.

Hi @charlesdarwin yes you are right its not a regular vcf, I should probably just refer to it as a tab-delimited file. No I am not sure that I would get that, I am not entirely sure what format indels are supposed to be in, but I do feel like its super weird for them to have more than one line associated with a single position, and eventually I need to change this into a VCF file and I feel like that would probably mess things up.

Did they perform exome sequencing or use a SNP chip? This is likely to occur for the latter, since they're then reporting every probe.

@DevonRyan I believe they used a SNP chip.

The VCF format does allow multiple variants at the same chromosome & position - in fact the official dbSNP VCFs have exactly that.

First of all, take these data with a very large pinch of salt. This sort of targeted analysis is not designed to produce high quality genetic data but to give an idea of a sample's ancestry. Given a file like the one you describe with no quality information, you have no way of knowing how reliable any of the variants called actually is. From personal experience with another, similar company, I expect that they aren't reliable at all.
I have tested four individuals using a commercial tool like the one you describe and while I don't doubt the validity of their main results (those displayed on their website), their raw data showed combinations of variants that were extremely unlikely to be real given that the individuals were still breathing.
This isn't to say that such companies don't do their job well, only that their job is limited to what they actually display. When you look into their raw data, you're on your own. Presumably, they apply various filters before displaying the main results, but these are most likely not applied to the raw data. So if you want to do your own analyses, this is unlikely to be a reliable dataset. Also, I would strongly urge you to seek the help of an expert when interpreting any data you do have.

That said, the rsIDs are simply identifiers of a specific genetic variation, not of a location. Any genomic position can have many (essentially infinite) possible variations. For example, a C residue at position 100 of chromosome 1 could become a G, a T or an A. Each of these will have its own, separate rsID. Now, if you also consider insertions and deletions, you can see why you can have essentially limitless possible variations at that specific position. So no, having more than one rsID describing the same position is not strange at all.
The way these tests work is they look at a list of possible variants, each of which has its own rsID, and they report what they find for that rsID in the sample. Then, they will report what they find at the rsID's position for each allele. This is explained in the FAQ article you linked to:

The possible observations are A for adenine, C for cytosine, G for guanine, T for thymine, I for insertion’ and ‘D for deletion, or 0 for missing data. Column one provides the identifier (including the #rsID where possible). Columns two and three contain the chromosome and basepair position of the variant using human reference 37.1 coordinates. Columns four and five contain the two alleles observed at this variant (genotype)

Let's have a look at the two rsIds you mention.

rs6010164: This is a single nucleotide variant (SNV) which describes a C residue where the reference genome has a T.

rs267607236: this is a deletion at position 50076841, so a very different variant to the above. In this one, the sequence TGAG is deleted.

Your results are:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

So, your sample doesn't have rs6010164 since it has the same sequence at that position as the reference. However, it does have the insertion as indicated by the two Is.
The Is and Ds of the insertions and deletions (indels) are going to be the easy bit here. The only way you can parse this sort of file and convert it to VCF is to get the list of target SNV rsIDs, check what they describe as the reference and as the variant and then, for every one, check which of the two is shown in your file.
If you're not an expert in this field, you will probably be better off hiring someone to do this for you. There are actually companies that offer this sort of service and can take these tab-delimited files as input. I even work for one!

It's pretty common for SNPs and indels that map to the same position in the genome to be assigned different rsIDs, particularly when the indel spans multiple base pairs in the reference. You can get more info about these particular variants here and here. Or you can see how both of these variants relate to the reference sequence here.
Without info about the reference allele it's pretty confusing (particularly because they appear to be reporting the opposite strand from the database (G vs. C/T), but I would interpret this to mean that both alleles match the reference for rs267607236 (they have the insertion), but do not match the reference for rs6010164 (G-C pair rather than T-A). This is also the most common genotype at this position, being found in almost 80% of people in most populations.

Thank you @heathobrien! If I want to convert my data to a VCF don't I need to one use only one of the rsid's since they are at the same position? Is there any standardized way to pick which one to use? And sorry, my bad, they did report from the forward strand.  I had changed the nucleotides since its actually not my own DNA data I'm posting here, but I should have checked to make sure I was changing it to a genotype that made sense.

You need to understand a few key points for your example given:
(a) Ancestry (and 23andMe) will use the letters I and D to indicate it is an InDel/DELINS. You have to read the rsID definition to determine if the variant is an Insert or Delete from the reference.  And therefore which letter is used to mean positive / variant and which for negative / reference.  They simply use the other letter to indicate negative / reference.
rs267607236 is defined as a deletion. So the I's mean the reference was detected in this case. Therefore the Ancestry file indicates no deletion of the 4 base-pairs at coordinates 50515271-4.
(b) Ancestry and dbSNP VCF annotation files uses a "left-aligned" coordinate for ALL InDel's and generally designate the coordinate BEFORE the variant.
This makes sense for Inserts. It gives the coordinate of the last base-pair in the reference just before the insert starts. That coordinate plus 1 is the next base-pair after the insert in the reference.
This same method is used for deletions.  But because the deleted values have coordinates in the reference, this makes deletions appear "one off". You could use the starting coordinate of the deletion when defining the variant (e.g. [TGAG/-]). Otherwise, you need the coordinate before if defining the deletion as [CTGAG/C].
rs267607236 is listed in the OP's information as coordinate 50515270. Because it is a deletion, the deleted sequence actually starts with the base-pair at coordinate 50515270+1 or 50515271.
(c) As it ends up, for rs267607236 in question here, the dbSNP VCF and the dbSNP online website entry are different (both latest release 156 in this case). This could happen but still be consistent because there is a T flanking the 4 base-pair deletion.  The starting reference sequence ...C TGAGT A... becomes CTA after the deletion.  But which T remains is not really known. The VCF file defines the deletion by retaining the trailing T. The online website retains the leading T.
And in fact, it appears the website entry starting coordinate is incorrect.  Given the way they defined this variant within the flanking sequence, the online website should use the coordinate 22:50515270 (like in the Ancestry file) and not 22:50515271 as it does now. This may have existed with that coordinate in an earlier release used by Illumina / Ancestry when defining the pipeline.  But why the VCF and online website are different now is confusing.
There are numerous other nuances, conventions and actual errors that have to be understand to properly convert a microarray raw data file to a VCF.
Here are all the entries that we referenced / consulted brought together for your perusal:
rs6010164      22    50515270    T    T
rs267607236    22    50515270    I    I

GRCh37 Reference sequence around rs267607236:
TTAC TGAGT AAGA  (70-74)

dbSNP Release 156 GRCh37 VCF:
NC_000022.10    50515270        rs6010164       T       C
NC_000022.10    50515269        rs267607236     CTGAG   C       # Note: leading T deleted

dbSNP web page definition:
Alleles:T>C             for rs6010164)
TTAC [T/C] GAGTAAGAT
22:50515270 (GRCh37)
NC_000022.10:50515270:T:C

Alleles:TGAGT>T         for rs267607236
TTACT [GAGT/-] AAGAT                                            # Note: trailing T deleted
22:50515271 (GRCh37)                # for consistency with flanking sequence, this should be ..270
NC_000022.10:50515271:TGAGT:T

It is not clear to me what Ancestry does when the same rsID has an insert and a deletion at the same position. For example, the line from dbSNP VCF for rs768019142 is:
NC_000001.11    10055   rs768019142     TA      T,TAA
which says there are two variants in the rsID. (This is common, especially for SNPs with multiple variants possible).  A deletion at 10056 OR a single letter insert after 10056. To represent this, instead of a single letter for the reference, they have given two base-pairs as the reference sequence. They probably should have defined two rsIDs at two adjacent locations. UCSC's versions of the VCF files are filtered and cleanup to simplify items like this.

For different variation type at the same location, each variation type is being assigned its own rsID by the dbSNP maintainers.
In your example, one chromosome in the pair has a single-nucleotide variation (T allele), while the other chromosome has an insertion variation.
Quote:

dbSNP only merges rs numbers that have an identical set of mappings to a genome and the same allele type (e.g. both must be the same variation type and share one allele in common).  We therefore would not merge a SNP and an indel (insertion/deletion) into a single rs number (different variation classes) since they represent two different types of mutational "events".

In your case, rs6010164 is of the Single Nucleotide Variation (SNV) variation type, while rs267607236 is of the Small deletions or insertions (DELINS) variation type.
Quote:

Despite its name, RefSNP is assigned to all variation types listed below with precise locations for both common and rare variations, including mutations. Most are typically small variations (<= 50bp).

Single nucleotide variation (SNV)
Short multi-nucleotide changes (MNV)
Small deletions or insertions
Small STR repeats
retrotransposable element insertions

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


text:
Gathering data on bacterial organism growth conditions
let's say I have a list of organism names like the example below
Achromobacter xylosoxidans
Acidithiobacillus ferrooxidans
Chloroflexus aurantiacus
Clostridium perfringens
Aquaspirillum arcticum

The goal is to identify the optimal growth conditions for each organism--meaning the optimal growth temperature and pressure. Obviously I can just google this information if the list is composed of a few organisms or even few tens of organisms. However, assuming that the list is too large, let's say several 100s of names, what would be an efficient means to gathering this information.
In searching for this info, I have not come across a database that has the information that I'm looking for in one place. Instead I find that this information is scattered across multiple sites, databases, scientific publications, protocols, etc. For each organism, I basically have to search through and extract this manually, and I was wondering if anyone can suggest some way to automate this. I am familiar with Python3, ideally I would accomplish this with python.
I though about implementing some sort of a web scraping solution, where I provide the name of the organism, as well as a few key words and see what happens. However, many times this data is not well tabulated and the name and whenever you do find information on the organism, you need to read to be able to interpret the information to get your answer...
maybe I'm not looking in the right place, or my approach to this problem is wrong/backwards/etc, any suggestions will help.
I haven't been able to find a centralized place to gather this information, that's why I posted here. The closest thing I've found is the ATCC.

https://www.biostars.org/p/235851/

yuo need to parse all of their product sheet, colored pdfs .....https://www.atcc.org/api/pdf/product-sheet?id=49176

You could use the BacDive API, or you could download the kgx formatted files from KG Microbe.

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


text:
UniProt flat file. Interpreting format of SUBCELLULAR LOCATION. What does flag mean?
In creating a parser for the UniProt flat file uniprot_sprot.dat (directory) in the manual for the section SUBCELLULAR LOCATION is listed

The format of SUBCELLULAR LOCATION is:
CC   -!- SUBCELLULAR LOCATION:(( Molecule:)?( Location.)+)?(
  Note=Free_text( Flag)?.)?
Where:

Molecule: Isoform, chain or peptide name
Location = Subcellular_location( Flag)?(; Topology( Flag)?)?(; Orientation( Flag)?)?

Subcellular_location: SL-line of subcell.txt ID-record
Topology: SL-line of subcell.txt IT-record
Orientation: SL-line of subcell.txt IO-record

What does Flag mean?
My best guess at present is that these are the text surrounded by parenthesis, e.g.
 Note=Not associated with microtubules (PubMed:21558460).

but that does not help in understanding why it is named flag.
While I can use my parser to verify this is correct when this section is working, it would be nicer to have a better understanding of what the word flag means in this context.
I also checked the XML schema but no help there.

Per request.
More examples from file.
These where gathered using grep.
grep -A 3 -n 'CC   -!- SUBCELLULAR LOCATION'  /mnt/.../uniprot_sprot.dat

71:CC   -!- SUBCELLULAR LOCATION: Host membrane {ECO:0000305}; Single-pass
72-CC       membrane protein {ECO:0000305}.

2163:CC   -!- SUBCELLULAR LOCATION: Membrane {ECO:0000305}; Single-pass membrane
2164-CC       protein {ECO:0000305}.

6003:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

7039:CC   -!- SUBCELLULAR LOCATION: Cell membrane {ECO:0000305}; Lipid-anchor,
7040-CC       GPI-anchor {ECO:0000305}. Note=In microneme/rhoptry complexes.
7041-CC       {ECO:0000250}.

7409:CC   -!- SUBCELLULAR LOCATION: Secreted {ECO:0000305}.

8677:CC   -!- SUBCELLULAR LOCATION: Virion {ECO:0000269|PubMed:9499098}. Host
8678-CC       endoplasmic reticulum-Golgi intermediate compartment
8679-CC       {ECO:0000269|PubMed:15016891, ECO:0000269|PubMed:9499098}.

12871:CC   -!- SUBCELLULAR LOCATION: Nucleus {ECO:0000250|UniProtKB:P48349}.
12872-CC       Cytoplasm {ECO:0000250|UniProtKB:P48349}. Note=Translocates from
12873-CC       the cytosol to the nucleus when phosphorylated.
12874-CC       {ECO:0000250|UniProtKB:P48349}.

grep -A 6 -n -m 200 'CC   -!- SUBCELLULAR LOCATION.*(.*)'  /mnt/.../uniprot_sprot.dat

1560049:CC   -!- SUBCELLULAR LOCATION: Secreted. Note=Dense core vesicles (DCV).

2452012:CC   -!- SUBCELLULAR LOCATION: DNA-(apurinic or apyrimidinic site) lyase,
2452013-CC       mitochondrial: Mitochondrion. Note=Translocation from the
2452014-CC       cytoplasm to the mitochondria is mediated by ROS signaling and
2452015-CC       cleavage mediated by granzyme A. Tom20-dependent translocated
2452016-CC       mitochondrial APEX1 level is significantly increased after
2452017-CC       genotoxic stress (By similarity). The cleaved APEX2 is only
2452018-CC       detected in mitochondria. {ECO:0000250}.

9243083:CC   -!- SUBCELLULAR LOCATION: Isoform Soluble copper chaperone CopA(Z):
9243084-CC       Cytoplasm {ECO:0000250|UniProtKB:Q59385}.

38197671:CC   -!- SUBCELLULAR LOCATION: E(rns) glycoprotein: Host cell membrane
38197672-CC       {ECO:0000269|PubMed:28290554}; Peripheral membrane protein. Virion
38197673-CC       membrane {ECO:0000269|PubMed:28290554}; Peripheral membrane
38197674-CC       protein {ECO:0000305}. Note=The C-terminus membrane anchor of Erns
38197675-CC       represents an amphipathic helix embedded in plane into the
38197676-CC       membrane.

4924705:CC   -!- SUBCELLULAR LOCATION: Note=(Microbial infection) In the presence
4924706-CC       of M.tuberculosis EsxA-EsxB complex decreased amounts of B2M are
4924707-CC       found on the cell surface (PubMed:25356553).
4924708-CC       {ECO:0000269|PubMed:25356553}.

02/18/2019
For the DE line it is noted

A block of DE lines may further contain multiple Includes: and/or Contains: sections and a separate field Flags: to indicate whether the protein sequence is a precursor or a fragment:

At present this seems to be related to the answer. When I have the file correctly parsed, then I can check for a correspondence as verification.

Of interest: UniProt subcell.txt - Controlled vocabulary of subcellular locations and membrane topologies and orientations.

Of interest: UniProt Evidence - Most information in UniProtKB has one or several “evidence tags” which describe the source of the information, e.g. an experiment that has been published in the scientific literature, an orthologous protein, a record from another database, etc.

Thanks. And sorry, I forgot to delete my comment after I added the example. Have you looked at the Swissknife perl module? This is produced by Uniprot precisely to deal with these files. That should at least point you in the right direction.

Ah, maybe. Perhaps to specify that the precursor is in one subcellular location while the mature product is in another. I don't know. And I couldn't find an example of these "flags" in an actual CC section of the uniprot_sprot.dat I had lying around. But I admit I just looked at a few, I don't really know how to grep for it.

Please don't update the question with what you learned, but post an answer instead! That way, the question will be marked as answered and can help future users. And thank you, I am very familiar with grep, I've been using it every day of my life for almost a couple of decades now.  :)

Thank you for asking this question. "Flag" is actually an obsolete notation which predates the introduction of proper evidence attribution in UniProtKB, and which escaped us in our efforts to update documentation to use the correct terminology. The correct  term would now be "Evidence", and the next version of the user manual will have this corrected (release 2019_06 due on July 3rd, 2019).
Before the current way of evidencing annotations was introduced in UniProtKB, the following "non-experimental qualifiers" or "flags" were used: By similarity, Potential, Probable. Absence of such a flag would generally imply that the information is based on an experiment.
When in doubt, please don't hesitate to contact the UniProt helpdesk.


text:
Are these standard species abbreviations and how to look up others?
Reactome uses a three letter code for each species which are listed below.
Are the letter codes used by Reactome a standard?
If so is there a web page or download file that can validate the code and translate it to a common name?
 HSA  Homo sapiens
 SSC  Sus scrofa
 GGA  Gallus gallus
 DME  Drosophila melanogaster
 RNO  Rattus norvegicus
 DDI  Dictyostelium discoideum
 MMU  Mus musculus
 ATH  Arabidopsis thaliana
 TGU  Taeniopygia guttata
 DRE  Danio rerio
 OSA  Oryza sativa
 BTA  Bos taurus
 XTR  Xenopus tropicalis
 CFA  Canis familiaris
 SCE  Saccharomyces cerevisiae
 PFA  Plasmodium falciparum
 CEL  Caenorhabditis elegans
 SPO  Schizosaccharomyces pombe
 ROT  Rotavirus
 FLU  Influenza A virus
 HCV  Hepatitis C Virus
 HBV  Hepatitis B virus
 HER  Human herpesvirus 1
 NME  Neisseria meningitidis serogroup B
 NME  Neisseria meningitidis
 CTR  Chlamydia trachomatis
 STY  Salmonella typhimurium
 STY  Salmonella enterica
 ECO  Escherichia coli
 OCU  Oryctolagus cuniculus
 XLA  Xenopus laevis
 CAE  Cercopithecus aethiops
 CGR  Cricetulus griseus
 HER  Human herpesvirus 8
 MCV  Molluscum contagiosum virus subtype 1
 MCV  Molluscum contagiosum virus
 HCY  Human cytomegalovirus
 HIV  Human immunodeficiency virus 1
 PCH  Penicillium chrysogenum
 TAE  Triticum aestivum
 VRA  Vigna radiata var. radiata
 VRA  Vigna radiata
 CFS  Crithidia fasciculata
 FCA  Felis catus
 AMA  Arenicola marina
 MGA  Meleagris gallopavo
 NGO  Neisseria gonorrhoeae
 SAU  Staphylococcus aureus
 MVI  Measles virus
 HAL  Human alphaherpesvirus 2

EDIT
Not a complete answer as this only covers a subset of all organisms but at present better than nothing.
Using KEGG Organisms
https://www.genome.jp/kegg-bin/get_htext
Enter letter code, e.g. GGA

EDIT - After reading answer by conchoecia
The graph database schema for Reactome only has two labels that use the 3 letter code which is named abbreviation in the schema:
Species and ReferenceTherapeutic
Species has taxid  (Taxon Identifier)
abbreviation     taxId       displayName
AMA              6344        Arenicola marina
ATH              3702        Arabidopsis thaliana
BAN              1392        Bacillus anthracis
BTA              9913        Bos taurus
CAE              9534        Cercopithecus aethiops
CAL              5476        Candida albicans
CBO              1491        Clostridium botulinum
CEL              6239        Caenorhabditis elegans
CFA              9615        Canis familiaris
CFS              5656        Crithidia fasciculata
CGR              10029       Cricetulus griseus
CPH              10703       Corynephage beta
CPO              10141       Cavia porcellus
CSB              60711       Chlorocebus sabaeus
CTE              1513        Clostridium tetani
CTR              813         Chlamydia trachomatis
CVI              10243       Cowpox virus
DDI              44689       Dictyostelium discoideum
DME              7227        Drosophila melanogaster
DRE              7955        Danio rerio
ECO              562         Escherichia coli
FCA              9685        Felis catus
FLU              11320       Influenza A virus
GGA              9031        Gallus gallus
HAL              10310       Human alphaherpesvirus 2
HAM              6706        Homarus americanus
HBV              10407       Hepatitis B virus
HCV              11103       Hepatitis C Virus
HCV              31649       Hepatitis C virus genotype 2a
HCV              31646       Hepatitis C virus subtype 1a
HCY              10359       Human cytomegalovirus
HER              10298       Human herpesvirus 1
HER              37296       Human herpesvirus 8
HGA              10376       Human gammaherpesvirus 4
HIV              11676       Human immunodeficiency virus 1
HPV              333760      Human papillomavirus type 16
HPV              337041      Alphapapillomavirus 9
HSA              9606        Homo sapiens
IBR              11120       Infectious bronchitis virus
LMO              1639        Listeria monocytogenes
LMO              1906951     Listeria monocytogenes serovar 1/2a
LPN              446         Legionella pneumophila
MCV              10280       Molluscum contagiosum virus subtype 1
MCV              10279       Molluscum contagiosum virus
MGA              9013        Meleagris gallopavo
MML              9544        Macaca mulatta
MMU              10090       Mus musculus
MTU              1773        Mycobacterium tuberculosis
MTU              83332       Mycobacterium tuberculosis H37Rv
MVI              11234       Measles virus
NGO              485         Neisseria gonorrhoeae
NME              491         Neisseria meningitidis serogroup B
NME              487         Neisseria meningitidis
OAR              9940        Ovis aries
OCU              9986        Oryctolagus cuniculus
OSA              4530        Oryza sativa
PCH              5076        Penicillium chrysogenum
PFA              5833        Plasmodium falciparum
RNO              10116       Rattus norvegicus
ROT              10912       Rotavirus
SAU              1280        Staphylococcus aureus
SCE              4932        Saccharomyces cerevisiae
SEV              11191       Sendai virus
SPO              4896        Schizosaccharomyces pombe
SSC              9823        Sus scrofa
STY              90371       Salmonella typhimurium
STY              28901       Salmonella enterica
TAE              4565        Triticum aestivum
TGO              5811        Toxoplasma gondii
TGU              59729       Taeniopygia guttata
VRA              3916        Vigna radiata var. radiata
VRA              157791      Vigna radiata
VVI              10245       Vaccinia virus
XLA              8355        Xenopus laevis
XTR              8364        Xenopus tropicalis

While ReferenceTherapeutic does not.
The species can be identified using the NCBI taxon identifier, e.g.
https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=6344

EDIT - After reading answer by user172818
Searching the file https://www.uniprot.org/docs/speclist.txt for 6344
AREMA E    6344: N=Arenicola marina
                 C=Lugworm
                 S=Lumbricus marinus

or using the URL: https://www.uniprot.org/taxonomy/6344

and with reference to

I think databases should stop inventing new species abbreviations.

The Wikipedia article on Lugworm shows the many Taxon identifiers in multiple databases.

No, there are no standardized three-letter abbreviations for species names. Three-letter abbreviations with roman letters is 26^3 combinations, which is only 17,576 possible combinations. As a result there are not enough three-letter abbreviations for all living species.
All three-letter abbreviations are either specific to a database or perhaps are shared between a database or two. There is currently no substitute for a binomial name (although NCBI taxid is a convenient parsable substitute for the species included in NCBI to date).

These are not. The closest to "standard" is the 5-character abbreviation by Swissprot. It names 25,886 species and has been used for decades. It is also easy to remember for some common species such as HUMAN, RAT, MOUSE and HORSE. I think databases should stop inventing new species abbreviations.


text:
Relationship between entries for same UniProt Id
I was playing around, looking at some AlphaFold predictions for certain proteins and I realized that sometimes the sequences of aminoacid look very different for the same UniProt id. For example, 15C8 from PDB has 2 chains:
DIVLTQSPAIMSASLGERVTMTCTASSSVSSSNLHWYQQKPGSSPKLWIYSTSNLA
SGVPARFSGSGSGTSYSLTISSMEAEDAATYYCHQYHRSPYTFGGGTKLEIKRADA
APTVSIFPPSSEQLTSGGASVVCFLNNFYPKDINVKWKIDGSERQNGVLNSWTDQD
SKDSTYSMSSTLTLTKDEYERHNSYTCEATHKTSTSPIVKSFNRN

and
EVQLQQSGAELVKPGASVKLSCTASGFNIKDTYMHWVKQKPEQGLEWIAQIDPANG
NTKYDPKFQGKATITADTSSNTAYLHLSSLTSEDSAVYYCAADPPYYGHGDYWGQG
TTLTVSSAKTTPPSVYPLAPGSAAQTNSMVTLGCLVKGYFPEPVTVTWNSGSLSSG
VHTFPAVLQSDLYTLSSSVTVPSSTWPSETVTCNVAHPASSTKVDKKIV

At the same time, when looking at the UniProt database (the corresponding id is P01869), there's a list of multiple structures and their 3D shapes (here: https://www.uniprot.org/uniprot/P01869). However, looking at the AlphaFold entry, I see that the amino acid sequence looks different, eg. for the case above it's:
AKTTPPSVYPLAPGSAAQTNSMVTLGCLVKGYFPEPVTVTWNSGSLSSGVHTFPAV
LQSDLYTLSSSVTVPSSPRPSETVTCNVAHPASSTKVDKKIVPRDCGCKPCICTVP
EVSSVFIFPPKPKDVLTITLTPKVTCVVVDISKDDPEVQFSWFVDDVEVHTAQTQP
REEQFNSTFRSVSELPIMHQDWLNGKEFKCRVNSAAFPAPIEKTISKTKGRPKAPQ
VYTIPPPKEQMAKDKVSLTCMITDFFPEDITVEWQWNGQPAENYKNTQPIMNTNGS
YFVYSKLNVQKSNWEAGNTFTCSVLHEGLHNHHTEKSLSHSPGLQLDETCAEAQDG
ELDGLWTTITIFISLFLLSVCYSAAVTLFKVKWIFSSVVELKQTLVPEYKNMIGQA
P

My question (as someone who's just getting started with proteins) would be: what is the relationship between all these entries corresponding to the same UniProt Id?

PDB structure 15C8 corresponds to the protein product of the mouse Igkc gene, with the UniProt accession P01837 as you can see in the linked page:

Your P01869 UniProt accession is for the protein product of the mouse Ighg1 gene instead, so why would you expect these to be the same sequence? They are different proteins coded for by different genes.

Thanks for your response! I got the P01869 after looking it up in this tool: https://www.uniprot.org/uploadlists/ (from PDB to UniProtKB). This returns both P01837 and P01869 -- do you happen to understand why?

And in any case, my confusion regarding the amino acid sequence still holds using the correct UniProt id you mentioned, as AlphaFold has the following sequence: RADAAPTVSIFPPSSEQLTSGGASVVCFLNNFYPKDINVKWKIDGSERQNGVLNSWTDQDSKDSTYSMSSTLTLTKDEYERHNSYTCEATHKTSTSPIVKSFNRNEC which is a substring of the first chain in my original post + "EC" at the end

@Andreea that's... a different question but yes, the structure predictions of AlphaFold or SwissFold or any other such collection tend to be of only parts of a protein and not the entire thing. As for why you got both, that's because P01837 and P01869 are similar proteins and the specific sub-sequence whose structure has been predicted occurs in both of them.


text:
Gene expression signature (AUC) from gene expression data with 48 +48 samples
Does anyone have a good method to find a signature where you combine the expression from multiple genes to predict a specific condition.
Assume you have performed sequencing of 48 cancer samples and 48 normal sample. Then find a combination of genes that predicts the cancer samples, using for instance AUC.

One approach is to use Elastic Net regression, I've seen lots of examples in the literature, though the methods a pretty thin in terms of implementation

I think that with cancer that is quite diverse and with this low number of samples you reproduce whatever gene signature in a larger cohort of samples. (Calculate the power of your study)

Rather than going with a specific (bioinformatic) package you could turn to more general machine learning approaches. For similar problems of mine (with much fewer samples), default implementations of Random Forests (which will also provide importance of features / genes) often yield perfect or near-perfect classification. Note that Random Forests have many favorable characteristics for working with genes, and also form the mechanistic basis of the best-performing (competition-winning) tools for the prediction of gene regulatory networks.

I have a method for doing this that uses a MCMC algorithm originally designed for population genetics (via the structure program) to generate a statistic between 0 and 1 indicating genetic risk. This needs the use of reference populations for cases and controls that the query individuals are added to for classification. I use T1D as an example. My methods paper is not peer-reviewed, and probably doesn't go into sufficient detail for what you are looking for, but at least indicates one way this can be done:
https://www.biorxiv.org/content/early/2017/01/31/104497
Based on the work I've done with this, using 48 cases and 48 controls with a very specifically defined type of cancer should be sufficient to create a generalised signature from human populations. If a generic cancer profile is needed, the populations are not representative of a general human population, or there are some cases in the control group, then things get trickier and require more individuals. On the other hand, 1500 per group is probably overkill.

There are many different ways of doing it! I would recomend the well established R pacakge CARET which have a whole chapter and many build in functions for it.
In addition you probably want to do some low-level filtering based on variance and fold changes (or the combination - ad DE analysis) to get the number of features down before progressing to the more advanced methods.


text:
How can the Autocovariances, autocorrelations, and autocorrelation coefficients be calculated from a Protein Amino Acid Sequence?
Given a normal protein sequence with the 20 standard amino acids, how can the 'Autocovariances', 'autocorrelations', and 'autocorrelation coefficients' of the sequence be calculated?
What is meant by these terms in the context of a protein sequence?
What are the appropriate algorithms and formulas?
Is there any existing package/function to calculate these values?

Could you provide some context? It could clarify the question and attract more answers. Maybe you read that in an article or it is for some specific goal in mind, or someone asked you this because a discussion about another topic

Protein autocorrelation is similar to usual autocorrelation of time series (that is, function $r(k)$ measuring correlation between values lagged by $k$), except that instead of time dimension we slide along the backbone of amino acid chain. The "dependent variable" can be any quantifiable property of the amino acids - then the autocorrelation measures similarity in that particular property.
See e.g. Fernandez et al., 2008, for more extensive background:
https://www.tandfonline.com/doi/full/10.1080/08927020802301920?scroll=top&needAccess=true


text:
knnImputation for missing data
I have a dataframe containing missing values. I want to impute them with knnImputation as shown below.
Now my question is: are data imputed row wise or column wise and is there a way to give this option in knnImputation?
A <-c(1,2,3,4,5, NA, NA)
B <-c(2,3,NA, 2,4, 7, 6)
C <-c(2,3,4, NA, 5, 1,7)
D<-21:27
data<-data.frame(A,B, C, D)
new_data<-knnImputation(data, k = 3, scale = T, meth = "weighAvg", distData = NULL)

According to the code it seems like by row scaling by column but I am not sure. But you can change the orientation of your matrix by transposing it with t . Another option is to test it with a symmetric (or quasi-symmetric) matrix

The data are imputed using values from the same column
Here's a simplified example to show what's going on:
A <-c(1,2,3,4,5)
B <-c(2,3,NA, 2,4)
C <- c(2,3,4, 6, 1)
data<-data.frame(A,B,C)
data

#>  A  B C
#>1 1  2 2
#>2 2  3 3
#>3 3 NA 4
#>4 4  2 6
#>5 5  4 1

# I'm using scale = F to make things simpler
knnImputation(data, k = 3, scale = F, meth = "weighAvg", distData = NULL)

#>  A        B C
#>1 1 2.000000 2
#>2 2 3.000000 3
#>3 3 2.594272 4
#>4 4 2.000000 6
#>5 5 4.000000 1

# calculate euclidean distances, dropping column of interest
dist(data[,-2])

#>1        2        3        4
#>2 1.414214
#>3 2.828427 1.414214
#>4 5.000000 3.605551 2.236068
#>5 4.123106 3.605551 3.605551 5.099020

# use k rows with smallest distance to row of interest to calculate weighted average
# using exp(-dist) as weights
weight_sum <- weight_sum = exp(-2.828427) + exp(-1.414214) + exp(-2.236068)
2*exp(-2.828427)/weight_sum + 3*exp(-1.414214)/weight_sum + 2*exp(-2.236068)/weight_sum

#>2.594271

If you want to impute using values from the same row, all you have to do is transpose the matrix:
t(knnImputation(t(data), k = 2, scale = F, meth = "weighAvg", distData = NULL))

#>     A        B C
#>[1,] 1 2.000000 2
#>[2,] 2 3.000000 3
#>[3,] 3 3.086729 4
#>[4,] 4 2.000000 6
#>[5,] 5 4.000000 1

(I had to lower k because there are only 2 extra columns)


text:
Reduce number of transcripts in a highly variable de novo transcriptome assembly
I have a de novo assembly using both multiple SRA and locally sequenced transcriptomes. I started with 270M PE reads from 9 tissues. Here are the assembly stats generated with TrinityStats.pl:
################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  543668
Total trinity transcripts:      1111089
Percent GC: 33.77

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2117
        Contig N20: 1324
        Contig N30: 945
        Contig N40: 714
        Contig N50: 561

        Median contig length: 344
        Average contig: 496.17
        Total assembled bases: 551293315

#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2172
        Contig N20: 1333
        Contig N30: 925
        Contig N40: 681
        Contig N50: 523

        Median contig length: 323
        Average contig: 473.63
        Total assembled bases: 257499729

As you can see, the total Trinity 'genes' number is very high. I suspect that this is due to a high biological variability of the transcriptomes (coming from organisms sampled in Italy, China, Spain...).

How can I assess whether the problem is really the biological variability?
How can I reduce this number? I'd rather not downsample the reads. It looks like CD-Hit might do what I'm looking for...

As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

Map RNA-seq reads to Trinity-generated transcripts using Salmon
Use the expression of BUSCO genes to set a credible signal cutoff
Subset transcripts based on this threshold [for us it was 50 counts]
Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:
cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

More details here:
https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9

Thank you for your answer. Prior to using CD-HIT, I decided to filter the reads according to BLAST hits with a custom pipeline. Only 150k transcripts had a hit, and of these, only 56K had a hit with the suitable taxon. This means that almost a million transcripts don't match with anything, and probably CD-HIT won't help much. Any idea on why this is happening? Also, given that BUSCO retrieved 975/978 core genes (so my assembly is basically complete), should I care about this number?

CD-HIT does an internal match; it doesn't care about matches to external databases. That BUSCO score is great; even C. elegans isn't 100%.

Yes, but my question is: should I worry about the fact that I have >1M transcripts, and of these only 5% seem to correspond to a mRNA? What's in the other 950k transcripts reconstructed by trinity?

Ah, and just for the record, today I ran CDHIT and it obtained 780k clusters, so yes, it improved the situation, but not to a great extent. This still looks unrealistic to me.

To better answer that question, I need more information about what you were sequencing, how the sequencing was done, and what sequences are coming up. For example, if there's a high chance of bacterial contamination and sequencing is for total RNA including microRNAs, then 780k is something that could happen.

Mind if we talk about it in chat? I'm afraid this is getting long for comments :)

Seems to me that your "genes" are very short. What sort organism (prokaryote, lower eukaryote, higher eukaryote) is this? I would look to see if your "genes" contained ORFs, and only keep those with sufficicntly long CDS. You will miss uORFs and ORF fragments this way, but you'd have much more confidence in what you did find.

It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values?

Higher eukaryote? Talking about giraffes? http://scienceblogs.com/evolgen/2006/10/18/do-higher-eukaryotes-smoke-mor/


text:
Correct way to assemble reference transcriptome - what's --samples_file option in Trinity?
I have 101 samples from 9 tissues of a nonmodel species, coming either from SRA or sequenced by my lab. I want to generate a reference transcriptome with Trinity for further analysis. In order to do this, I downsampled the reads in order to have 20 to 30M reads for each tissue, for a total of 270M PE reads, which I combined in two files, one for left reads and one for right ones.
Now,  I came across the --samples_file flag in Trinity and I'm not sure if I should use it instead of concatenating the transcriptomes. Any help on good practice on this?

--samples_file flag is for DE analysis, you can put info about biological replicates in that file. You say that you want to make a ref transcriptome, so you don't need a samples file for DE analysis.
Why do you down sample? Are you not afraid you will not have enough reads to have a complete transcriptome now?

Thanks, it's what I was thinking. // I have 270 million reads, I think it's enough. Using all transcriptomes means having >750 million reads, I don't have the computational resources to assemble it.

Okay, I understood 270 PE reads...

Yes, I actually wrote 270 PE reads, my bad.


text:
How to format the Nexus file to create geographical maps with haplotypes in PopArt?
I'm trying to create a map in PopArt that allows me to visualise the distribution of several species. These are my Traits (=species) and GeoTags blocks:
BEGIN TRAITS;
    Dimensions NTRAITS=9;
    Format labels=no missing=? separator=Comma;
    TraitLabels SNE Lsp A B C D 3n 4n Lep;
Matrix
SNE1 1,0,0,0,0,0,0,0,0
SNE2 1,0,0,0,0,0,0,0,0
SNE3 1,0,0,0,0,0,0,0,0
SNE4 1,0,0,0,0,0,0,0,0
SNE5 1,0,0,0,0,0,0,0,0
SNE6 1,0,0,0,0,0,0,0,0
PUR_4n 0,0,0,0,0,0,0,1,0
GRA 0,1,0,0,0,0,0,0,0
OJE 0,1,0,0,0,0,0,0,0
ESC 0,0,0,0,0,0,0,0,1
LNM 0,0,0,0,0,0,0,0,1
PIN 0,0,0,0,0,0,0,0,1
URD 0,0,0,0,0,0,0,0,1
VIS 0,0,0,0,0,0,0,0,1
LYE 0,0,0,0,0,0,0,0,1
BEN_A 0,0,1,0,0,0,0,0,0
ALJ1_A 0,0,1,0,0,0,0,0,0
ALJ2_A 0,0,1,0,0,0,0,0,0
ESP2_B 0,0,0,1,0,0,0,0,0
ESP3_B 0,0,0,1,0,0,0,0,0
ALC_C 0,0,0,0,1,0,0,0,0
AVA1_3n 0,0,0,0,0,0,1,0,0
AVA2_3n 0,0,0,0,0,0,1,0,0
PER_3n 0,0,0,0,0,0,1,0,0
TRI_3n 0,0,0,0,0,0,1,0,0
VDI_3n 0,0,0,0,0,0,1,0,0
VIL_3n 0,0,0,0,0,0,1,0,0
CVE1_4n 0,0,0,0,0,0,0,1,0
CVE2_4n 0,0,0,0,0,0,0,1,0
LNA1_4n 0,0,0,0,0,0,0,1,0
LNA2_4n 0,0,0,0,0,0,0,1,0
PRA1_4n 0,0,0,0,0,0,0,1,0
PRA2_4n 0,0,0,0,0,0,0,1,0
VEN1_4n 0,0,0,0,0,0,0,1,0
VEN2_4n 0,0,0,0,0,0,0,1,0
ZAR_4n 0,0,0,0,0,0,0,1,0
SFI1 0,0,0,0,0,0,0,0,1
SFI2 0,0,0,0,0,0,0,0,1
SIB 0,0,0,0,0,0,0,0,1
TIS1_D 0,0,0,0,0,1,0,0,0
TIS2_D 0,0,0,0,0,1,0,0,0
TIS3_D 0,0,0,0,0,1,0,0,0
TIS4_D 0,0,0,0,0,1,0,0,0
TVI_D 0,0,0,0,0,1,0,0,0
PVRFAS0001 1,0,0,0,0,0,0,0,0
PVRFAS0002 1,0,0,0,0,0,0,0,0
PVRFAS0024 1,0,0,0,0,0,0,0,0
PVRFAS0025 1,0,0,0,0,0,0,0,0
PVRFAS0027 1,0,0,0,0,0,0,0,0
PVRFAS0028 1,0,0,0,0,0,0,0,0
;
END;
Begin GeoTags;
    Dimensions NClusts=29;
    Format labels=no separator=Spaces;
    [ClustLatitude 37.15 37.77 37.13 39.37 37.11 40.60 40.59 37.04 40.61 38.71 36.70 43.48 40.18 40.24 40.57 40.26 37.10 37.86 38.52 36.56 40.23 38.01 39.58 38.75 40.54 40.55 38.70 40.04 40.69;
    ClustLongitude -3.70 -3.02 -3.52 -3.72 -3.03 -4.33 -4.33 -3.43 -2.44 0.05 -5.35 3.55 -6.55 -2.06 -4.16 0.26 -3.26 -1.56 -3.56 -4.86 -2.07 -2.87 -3.87 0.18 -4.18 -4.19 -0.49 -7.49 -2.59;]
Matrix
SNE1 37.044888 -3.135564
SNE2 37.044888 -3.135564
SNE3 37.044888 -3.135564
SNE4 37.044888 -3.135564
SNE5 37.044888 -3.135564
SNE6 37.044888 -3.135564
PUR_4n 37.133158 -3.524806
GRA 36.708272 -5.348794
OJE 36.562847 -4.853522
ESC 40.5659132 -4.2236186
LNM 40.609845 -4.344144
PIN 38.752157 0.180615
URD 39.370381 -3.720659
VIS 38.520059 -3.564201
LYE 39.580131 -3.870267
BEN_A 38.712603 0.052372
ALJ1_A 38.704010 -0.490801
ALJ2_A 38.704010 -0.490801
ESP2_B 37.865615 -1.569499
ESP3_B 37.865615 -1.569499
ALC_C 40.264366 0.267744
AVA1_3n 43.490194 3.563363
AVA2_3n 43.490194 3.563363
PER_3n 40.610461 -2.447806
TRI_3n 40.697087 -2.595488
VDI_3n 40.244599 -2.067766
VIL_3n 40.231843 -2.078962
CVE1_4n 40.559003 -4.199853
CVE2_4n 40.559003 -4.199853
LNA1_4n 40.597289 -4.330246
LNA2_4n 40.597289 -4.330246
PRA1_4n 37.113314 -3.035071
PRA2_4n 37.113314 -3.035071
VEN1_4n 37.155124 -3.703330
VEN2_4n 37.155124 -3.703330
ZAR_4n 40.547662 -4.184395
SFI1 40.042132 -7.498515
SFI2 40.042132 -7.498515
SIB 40.182598 -6.551259
TIS1_D 37.771593 -3.023286
TIS2_D 37.771593 -3.023286
TIS3_D 37.771593 -3.023286
TIS4_D 37.771593 -3.023286
TVI_D 38.0140003 -2.8751617
PVRFAS0001 37.100000 -3.260000
PVRFAS0002 37.100000 -3.260000
PVRFAS0024 37.100000 -3.260000
PVRFAS0025 37.100000 -3.260000
PVRFAS0027 37.100000 -3.260000
PVRFAS0028 37.100000 -3.260000
;
End;

If I upload my Nexus file like this, the program throws an error saying "Too many rows in GeoTags block." However, I have 50 samples in both Traits and GeoTags. The same happens if I uncomment (in Nexus format, comments are in brackets) the ClustLatitude and ClustLongitude files.
If I set the "labels" flag on "yes", it takes the traits as the geographical locations and tries to find locations in the world starting with the trait names.
If I comment the "GeoTags" block, PopArt correctly outputs a network.

What am I missing? Does somebody have a Nexus example where the traits do not correspond to geographical location of the samples?

I don't know this format but you seem to have two BEGIN and two ENDs. To me it looks like you have two datasets in one. Also one is in caps and the other is not. Are you sure that the file is in the right format? Did you check it (if there is any tool to do that) ?

This program uses a custom version of the Nexus format. Here is the manual: http://popart.otago.ac.nz/doc/popart.pdf and it uses exactly the "begin" and "end" blocks I'm using - weird format, I know, but I'm trying to stick as much as possible to the example given...

@Llopis Additional info here (Supplementary data of the PopArt paper) and here (Nexus example in the PopArt site)

In the pdf you link in the first comment the traits data is formatted different (see the corresponding section of "Traits Tables"). Are you sure you should provide all these data in a single file? The error seems to imply that when an END is found it should be the end of the file or something alike.

Ah you mean the indent? That's a copy-paste error, I'll edit that. In the original file I indented the block between "BEGIN TRAITS" and "Matrix". Note that if I comment the GeoTags block, a network is produced correctly. The manual notes: " I've separated the individual blocks with text describing the important points, but they should all be included in a single text file, as shown in here."

No, I mean that the section "traits Tables" describe a different format (coma separated values I think) than the one you have, but trait information is described in two formats.

authors: Обедающий философ


text:
Merging transcriptomes coming from different experiments
I'm planning to build a transcriptome by pooling all existing transcriptomes in SRA for a non-model species (which has no reference genome) to study differentially expressed genes and the like. It might worth mentioning that the transcriptomes I want to pool mostly come from different tissues and sampling sites. These transcriptomes would come from five SRA BioProjects plus one project from my laboratoryfor a total of 102 single transcriptomes coming from about 20 individuals (plus three pooled larvae transcriptomes). I'm especially interested in differentially expressed genes in different tissues, on genes specific to different tissues and on the UTRs of differentially expressed genes.
My idea is to build a common reference transcriptome against which I'd re-assemble the individual transcriptomes, and only then I'd get differential expression data etc.
I'm looking for literature (or personal experience) dealing with assemblies of transcriptomes / comparing transcriptomes coming from different datasets in order to understand whether this is a good idea or not and get a grasp of what technical and biological problems I could face. I'm (vaguely) aware of the batch effect and of one way to remove it, but I'd like to read papers dealing with my same problem (different SRA databases, without reference genomes). Up until now I found just projects dealing with microarray data.

R's Bioconductor (obviously). There is a lot of experience on this board in exactly this area using NGS, albeit with the focus on "model" species (humans). No reference genome, sounds a bit scary.

I don't have personal experience, so I won't post an answer. Be aware the with ComBat you could introduce a new batch effect (it is not necessary to remove them in order to work with them, you can take it into account). The other problem is that you will have several conditions (even the same tissue is not the same at the first year of the organism that just a year before dying).

However the most important thing is what do you expect to learn? The differentially expressed genes don't mean anything if you don't know what you are comparing or if you are not looking for something very specific (healthy vs disease, one tissue vs the other...). Please [edit] the question to explain what is your biological question and it will be easier to help. (Also if you provide the number of projects of your organism and the number of total samples it would be helpful to us too)

@llrs I edited the question with the things you asked and I clarified it a little bit (I was too much into my project to notice that I had skipped some fundamental information...). I'm aware I will have different conditions, that's why I'm not entirely sure that this is a good idea (i.e. that can give me some significant biological insight).

@LinuxBlanket Ok, that information focused the question. Still you haven't explained what are you trying to do with this comparison. Why do you want to compare transcriptomes coming from different datasets? Having several conditions is not a problem if you want to identify what is more important of each condition. But I'm not sure what comparison are you willing to do

@llrs I'm especially interested in genes expressed exclusively in gonads, but I also want to have some insights about the basic transcriptomic features of this animal, like its housekeeping genes. That's why I need to have as much tissues as possible, so that I can really know which mRNAs are transcribed only by the gonads and which ones are transcribed by every tissue instead.

Great! (Do not forget to [edit] this into the question.) How many samples do you have from the gonads and from other tissues? If this is your objective and this is the best you can find I do not see any problem in doing it. The biological problems are hard without knowing the organism. The technical problems can't be addressed without finding them... so I'm not sure if I can answer more

As mentioned in the comments, your biological question is going to be vitally important in determining the best way to tackle this issue. However, there are some things that you should generally keep in mind.

Most de novo assemblers will perform kmer normalization. This means that the merged RNA-Seq data will be normalized across all libraries. This is almost certainly not what you want because you'll be randomly selecting reads from all libraries without regard to their treatments. So, if you're going to allow for kmer normalization, you should normalize prior to merging the libraries so that you have equal representation of the reads by library. If you don't normalize, you're going to have a very fragmented assembly, and it will take a long time to generate.
De novo assemblers using short reads assume all identified exons with junction support may be used to generate splice variants with equal probability. I'm not entirely sure how to best phrase that, but here's an example. Suppose in a cold-treated library that Gene A generates a transcript using exons 1-2-3, but a heat-treated library generates transcripts using 3-4-5 and 2-3-5. It may be that merging the libraries results in transcripts generated containing various permutations of those exons (e.g. 1-2-3-5, 2-3-4, etc.) that aren't generated in any of the biological settings but can no longer be distinguished. This will lead to false positive transcripts which will detrimentally affect expression quantification.
Merging assemblies is difficult. So, to be true to each treatment, it makes more sense to generate individual assemblies (which will also result in library-specific kmer normalization). StringTie essentially does this using reference genomes. You assemble using each library (all bio-reps for a treatment), and then merge after the fact and re-quantify. This is the best approach. However, it can be difficult and computationally taxing. Some programs like CD-HIT (specifically cd-hit-est) are typically used in these cases. Keep in mind that CD-HIT will drop sequences with high similarity, not merge them. This may lead to loss of some uniquely assembled sequence. Other approaches involve using programs like CAP3 (a DNA assembler) to treat the transcripts like longer reads for assembly. Depending on the parameter choice, this might work, but you will risk generating chimeric transcripts.

How you move forward with the assembly process will again depend on your biological questions. I would tend to promote that you tackle the questions one at a time rather than using a single transcriptome to address them all, but there may be circumstances in which merging assemblies makes more sense.
In your specific case, I would perform multiple assemblies using a variety of the approaches described above and then annotate the assemblies (ex. using Trinotate) to see how their quality metrics compare. I would discourage generating a single transcriptome and then re-assembling individual transcriptomes from that assembly -- the single transcriptome may have misassemblies that would then propagate to the individual assemblies, and this process would generally be messy/difficult.

Thank you for your answer. I edited my question to clarify it a little bit: I want to use the assembly as a reference to map reads from the single transcriptomes and do (among others) DE analysis. I probably won't use the sequence information directly from the reference.

As mentioned above, "The biological problems are hard without knowing the organism. The technical problems can't be addressed without finding them." I would second that statement and direct you to the last paragraph I just added for clarification.


text:
Why do reverse radtags have different start points in radtag sequencing?
Why is it that when doing paired end Radseq that R1 and R2 have different start points? For instance, while doing Radseq you do a restriction digestion of your DNA using enzyme X- shouldn't the sequencing reaction take place from either end of where the enzyme cut? If it doesn't- but happens randomly where do your multiplex sequences go?
For instance- if you are sequencing the sequence on an illumina machine with the "<>" representing the cut site-
Where does the sequencing reaction start from in each direction?
Where does the multiplexing sequence go? - does it go on each?
5-GGCTAGC<CCTGCA*GG>GCGCTT-3
 3-CCGATCG<GG*ACGTCC>CGCGAA-5

Hi @jawz_34, it appears you are asking several questions, and some of the questions are not relevant to the question title. It will be much easier for others to answer your post if you structure your post such that all of the post is remains relevant to the post title.

Assume that after digestion you're left with the following sequence:
5' AGCTAGCTAGCTAGCTGACGTCGTAC 3'
3' GTACGACGTCAGCTAGCTAGCTAGCT 5'

Then one read would start with AGCTAG and the other would start with TCGATC.
The barcodes (what are used for multiplexing) are outside of this and part of the illumina adapters. They're sequences in completely different reads. This is nicely illustrated below:

Index 1 and 2 don't actually need to be present. Even if they are, whether the sequencer actually sequences them is a user setting.
For those curious, the order of the reads as they're sequenced are:

read #1
Index 1 (optional)
Index 2 (optional)
read #2 (optional)

If you ever get more than 2 reads from a sequencing provider, then it's likely that they included the index sequence (I have no clue why, but some sequencing centers like to do this...even though it confuses the heck out of their users).


text:
AND or OR gene - protein logical relation in KEGG
How do I know if it is a AND or OR relationship between genes associated with the same enzyme in KEGG? E.g., for EC:1.6.1.2 (http://www.genome.jp/dbget-bin/www_bget?ec:1.6.1.2) does the bacterium Streptomyces coelicolor (SCO) need SCO7622 AND SCO7623 or SCO7622 OR SCO7623?

The enzyme EC:1.6.1.2 is an NAD(P) transhydrogenase and consists of two subunits. The two subunits are encoded by two different genes pntA and pntB.
So in terms of computer logic I would say AND, both subunits are necessary for a functional NAD(P) transhydrogenase enzyme (EC:1.6.1.2).

This was only an example, is there a way of knowing in general?

With enzymes from KEGG there is not, I am afraid. In your example it becomes clear that we see 2 subunits, and hence form the same protein. However, if you look for example at human isocitrate dehydrogenase (IDH) EC 1.1.1.42, you will see that two different genes (IDH1 and IDH2) are listed. These 2 genes are not subunits of the same protein but are 2 different proteins (one in mitochondrium and the other in cytosol). With KEGG the interest in enzyme, not gene.


text:
SPAdes error during assembly
I try to perform a hybrid assembly using Unicycler. Unicycler use SPAdes to assembly the illumina sequences, and the program crashes while in SPAdes. The spades.log file says something about an invalid argument and error code 22. I have run the spades.py --test, and everything seems to work fine.
Command line: /home/user/SPAdes-3.11.1-Linux/bin/spades.py    -1  /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz -2  /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz -o  /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction --threads   8   --only-error-correction

System information:
  SPAdes version: 3.11.1
  Python version: 2.7.12
  OS: Linux-4.10.0-40-generic-x86_64-with-Ubuntu-16.04-xenial

Output dir: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction
Mode: ONLY read error correction (without assembling)
Debug mode is turned OFF

Dataset parameters:
  Multi-cell mode (you should set '--sc' flag if input data was obtained with MDA (single-cell) technology or --meta flag if processing metagenomic dataset)
  Reads:
    Library number: 1, library type: paired-end
      orientation: fr
      left reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz']
      right reads: ['/home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz']
      interlaced reads: not specified
      single reads: not specified
Read error correction parameters:
  Iterations: 1
  PHRED offset will be auto-detected
  Corrected reads will be compressed (with gzip)
Other parameters:
  Dir for temp files: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/tmp
  Threads: 8
  Memory limit (in Gb): 87

======= SPAdes pipeline started. Log can be found here: /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/spades.log

===== Read error correction started.

== Running read error correction tool: /home/user/SPAdes-3.11.1-Linux/bin/hammer /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info

  0:00:00.000     4M / 4M    INFO    General                 (main.cpp                  :  83)   Starting BayesHammer, built from N/A, git revision N/A
  0:00:00.001     4M / 4M    INFO    General                 (main.cpp                  :  84)   Loading config from /home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info
  0:00:00.045     4M / 4M    INFO    General                 (memory_limit.hpp          :  51)   Memory limit set to 87 Gb
  0:00:00.045     4M / 4M    INFO    General                 (main.cpp                  :  93)   Trying to determine PHRED offset
  0:00:00.071     4M / 4M    INFO    General                 (main.cpp                  :  99)   Determined value is 33
  0:00:00.071     4M / 4M    INFO    General                 (hammer_tools.cpp          :  36)   Hamming graph threshold tau=1, k=21, subkmer positions = [ 0 10 ]
  0:00:00.072     4M / 4M    INFO    General                 (main.cpp                  : 120)   Size of aux. kmer data 24 bytes
     === ITERATION 0 begins ===
  0:00:00.100     4M / 4M    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 289)   Building kmer index
  0:00:00.101     4M / 4M    INFO    General                 (kmer_index_builder.hpp    : 106)   Splitting kmer instances into 128 buckets using 8 threads. This might take a while.
  0:00:00.101     4M / 4M    INFO    General                 (file_limit.hpp            :  32)   Open file limit set to 1024
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  74)   Memory available for splitting buffers: 3.62484 Gb
  0:00:00.101     4M / 4M    INFO    General                 (kmer_splitters.hpp        :  82)   Using cell size of 524288
  0:00:00.102     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R1_001.fastq.gz
  0:00:57.279     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 1338594 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 2152556 reads
  0:03:30.505     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             :  96)   Processing /home/user/minion/organism/Illumina/BD_S3_L001_R2_001.fastq.gz
  0:09:23.290     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 3442412 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 106)   Processed 4305112 reads
  0:14:08.692     8G / 8G    INFO   K-mer Splitting          (kmer_data.cpp             : 111)   Total 4305112 reads processed
  0:14:08.726    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 109)   Starting k-mer counting.
  0:20:36.408    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 115)   K-mer counting done. There are 205776690 kmers in total.
  0:20:36.409    32M / 8G    INFO    General                 (kmer_index_builder.hpp    : 121)   Merging temporary buckets.
  0:24:43.505    32M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 298)   Building perfect hash indices
  0:25:40.125   100M / 8G    INFO    General                 (kmer_index_builder.hpp    : 137)   Merging final buckets.
  1:11:46.553   100M / 8G    INFO   K-mer Index Building     (kmer_index_builder.hpp    : 320)   Index built. Total 95426752 bytes occupied (3.70991 bits per kmer).
  1:11:46.598   100M / 8G    INFO   K-mer Counting           (kmer_data.cpp             : 357)   Arranging kmers in hash map order
  1:12:55.597     3G / 8G    INFO    General                 (main.cpp                  : 155)   Clustering Hamming graph.
  1:21:35.596     3G / 8G    INFO    General                 (main.cpp                  : 162)   Extracting clusters
=== Stack Trace ===
[0x407a6a]
[0x409561]
[0x40c79b]
[0x40129a]
[0x5442b0]
[0x40620d]
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
Verification of expression '(intptr_t) MappedRegion != -1L' failed in function 'void MMappedWriter::reserve(size_t)'. In file '/spades/src/common/io/kmers/mmapped_writer.hpp' on line 94. Message 'mmap(2) failed. Reason: Invalid argument. Error code: 22'.
hammer: /spades/src/common/io/kmers/mmapped_writer.hpp:94: void MMappedWriter::reserve(size_t): Assertion `(intptr_t) MappedRegion != -1L' failed.

== Error ==  system call for: "['/home/user/SPAdes-3.11.1-Linux/bin/hammer', '/home/user/minion/organism/hybrid_assembly/spades_assembly/read_correction/corrected/configs/config.info']" finished abnormally, err code: -6

In case you have troubles running SPAdes, you can write to spades.support@cab.spbu.ru
or report an issue on our GitHub repository github.com/ablab/spades
Please provide us with params.txt and spades.log files from the output directory.

Answer from @mjoppich, converted from comment:
Late answer, but it looks like you are running out of memory, hence the intptr is invalid.


text:
Finding out proteins that have a similar domain to a specific toxin
I am trying to find out that which DNA sequences(in the below) have a similar domain to a specific toxin which is Clostridium botulinum
You can find the fragment of the DNA sequence that belong to that toxin gene below:

gi|27867582 (fragment of the known Clostridium botuninum toxin gene)
  GTGAATCAGCACCTGGACTTTCAGATGAAAAATTAAATTTAACTATCCAAAATGATGCTT
  ATATACCAAAATATGATTCTAATGGAACAAGTGATATAGAACAACATGATGTTAATGAAC
  TTAATGTATTTTTCTATTTAGATGCACAGAAAGTGCCCGAAGGTGAAAATAATGTCAATC
  TCACCTCTTCAATTGATACAGCATTATTAGAACAACCTAAAATATATACATTTTTTTCAT
  CAGAATTTATTAATAATGTCAATAAACCTGTGCAAGCAGC

I ran a blastx job over the above sequence which translates nucleotide into protein. Here's the result.
The result tells me that. The particular DNA sequence listed above, belongs to botulinum neurotoxin type E [Clostridium botulinum]
Besides, I am given two DNA sequences:
>Sequence A
TCTATCAAGTAGATTATTAAATACTACTGCACAAAATGATTCTTACGTTCCAAAATATGA
TTCTAATGGTACAAGTGAAATAAAGGAATATACTGTTGATAAACTAAATGTATTTTTCTA
TTTATATGCACAAAAAGCTCCTGAAGGTGAAAGTGCTATAAGTTTAACTTCTTCAGTTAA
TACAGCATTATTAGATGCATCTAAAGTTTATACGTTTTTTTCTTCAGATTTTATTAATAC

> Sequence B
TCCTGGCTCAGGACGAACGCTGGCGGCGTGTGCTTAACACATGCAAGTCGAGCGATGAAG
CTTCCTTCGGGAAGTGGATTAGCGGCGGACGGGTGAGTAACACGTGGGTAACCTGCCTCA
AAGTGGGGGATAGCCTTCCGAAAGGAAGATTAATACCGCATAACATAAGAGAATCGCATG
ATTTTCTTATCAAAGATTTATTGCTTTGAGATGGACCCGCGGCGCATTAGCTAGTTGGTA

Again, I translated them into proteins. The result are listed respectively(A and B):
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z423E71C015
https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=Z42AKVWC015
The first job tells that, Sequence A belongs to botulinum neurotoxin type F [Clostridium baratii] which means that similar domain with the above botulinum neurotoxin>
Second job tells me that SequenceB belongs to unknown protein [Clostridium botulinum NCTC 2916]. So that means it is not in the same domain with the above botulinum neurotoxin right?

Looks like you've got two identical links for your blast searches

Ups. I edited the question.

What's a 'domain' here? Are you talking about protein domains? Domains of life? As in animals, plants, bacteria etc? What proteins were you given? You show DNA sequences here, not proteins. Please [edit] your question and clarify.

A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.
This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast results do include conserved domain information though.

The sequence B has no conserved domain, but sequenceA and first sequence (toxin gene) has one. When I click to the conserved domain information, they show the same result. That means they are in the similar domain, right?

actually, when I click "search for similar protein architectures" in Sequence A, I see the botulinum neurotoxin type E which is the protein in the original sequence.

Sounds like you have your answer


text:
Range of values for VAF in heterozygous genotype
(VAF = variant allele frequency)
This question follows from this one.
What is the expected range of values for VAF of mutations with heterozygous genotype (i.e., 0/1 or 0|1)?
Is it possible to have 0% and 100% (respectively having allele-specific depth of 0 in ALT or REF, with non-zero depth in the other allele)?

From the link that I edited into your question (realized that the answer was there):

VAF is the percentage of sequence reads observed matching a specific DNA variant divided by the overall coverage at that locus. Because NGS provides a near random sample, VAF is thus a surrogate measure of the proportion of DNA molecules in the original specimen carrying the variant. For constitutional genetic testing, VAF is a measure of diploid zygosity: heterozygous loci should be near 50% VAF, homozygous loci should be near 100%, and reference loci should be near 0%. Deviations from these three expected values should be considered suspicious as potential errors due to incorrect base calls or alignment. For somatic testing, contamination from normal cells and tumor heterogeneity combine to cause unpredictable VAFs. In some clinical scenarios, such as testing a patient for therapeutic resistance mutations, the desired sensitivity dictates that variants with low VAF be included in downstream analysis. This all means that in order to maintain an acceptable level of sensitivity, variant filtering for somatic variant must be highly permissive from a technical validity perspective.

This is assuming diploid humans. For other ploidies the expected proportions of heterozygosity will match the possible heterozygous genotypes (e.g. 33%/67% for triploid heterozygotes, 25%/75% for tetraploid, etc.).


text:
How to de novo hybrid assemble with Pacbio CCS and Illumina PE reads
I would like to perform de novo genome assembly on a diploid microalgal strain.
I have two datasets:

PacBio CCS/HiFi reads, low coverage.
Illumina PE 2x150 (standard shotgun)

Does anybody have any software recommendations for how best to use both in a hybrid de novo assembly?
Most options out there are set up for Pacbio CLR reads, not Pacbio CCS/HiFi, for using in hybrid assembly. Perhaps I should assemble short read contigs first then combine with wengan, or go back to the PacBio subreads instead of CCS to use in a hybrid method?

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

If you have reasonable coverage with HiFi (30X per haplotype), then there usually isn't any reason to also use Illumina. My understanding is that Illumina can make the assembly worse if you try to polish with it, because it can erroneously induce haplotype switching errors.
I would use HiCanu or HiFiAsm, depending on whether you expect significant heterozygosity, and also assuming that you have non-trivial HiFi coverage.
If your HiFi coverage is very low (<15X per haplotype, maybe), you might consider going back to subreads, assembling those, and polishing with Pilon and your Illumina reads, or one of the workflows that you have found. But otherwise I don't think the Illumina will be very useful. If you had Hi-C reads you could use them directly in HiFiAsm to assist phasing and scaffolding, but I'm assuming that your Illumina is shotgun.
I know it's a PB promo material, but you see that PB promotes HiFi as not needing polishing, as HiFi has accuracy >99%, so shotgun Illumina is not going to do much. There is some argument for still polishing with your original long reads.
I can recommend any of the linked papers for ideas or further explanation.

Unfortunately, my HiFi coverage isn't that great, (unexpectedly large genome). However, my shotgun short read dataset is large and good quality, hence the question.

@bishopia Hmmmm I see what you mean. In that case, subread assembly+Illumina polishing might be your best bet. It would still be helpful to know more specifics (e.g. coverage, read length profile of each dataset: CCS, subread, ILMN). I am pretty skeptical of the usefulness of Illumina contigs as the basis of assemblies of large/complex/heterozygous genomes, compared to even low-coverage PB data.


text:
Understanding regularization step in Meta-cell pre processing
I'm reading through this paper: Baran, Y., Bercovich, A., Sebe-Pedros, A., Lubling, Y., Giladi, A., Chomsky, E., ... & Tanay, A. (2019). MetaCell: analysis of single-cell RNA-seq data using K-nn graph partitions. Genome biology, 20(1), 1-19.

Maybe I am not familiar with this notation but I don't understand the regularization procedure as follows:

What is happening in these max functions? What does the 0 indicate? They have some terms $x$ and then in the max function they write $max( x, 0)$.
I assume $K$ is some regularization parameter, like a cost. Is $s_{ij} * s_{ji}$ the rank product?  I.e. the geometric mean?

Without looking in too much detail:

We may suppose that $S$ is a matrix with dimensions indexed as $(i, j)$. So $s_{ij}$ is the cell of $S$ indexed to the $i$th row and $j$th column, and $s_{ji}$ is the cell vice versa. Therefore, symmetrizing $S$ is a matter of ensuring that $s_{ij} == s_{ji}$, which you can achieve by multiplying the (ranked?) entries. Symmetric matrices are nice for doing certain kinds of math, and yield a handy representation of undirected graphs, which the authors appear to be using for kNN.
$max(x, 0)$ is just saying that you want either $x$ if $x$ is non-negative, or $0$ if it is negative. Thus, you make the whole matrix non-negative, which makes some handy matrix math possible. You can see how this works computationally by inspecting e.g. the Python max() function.


text:
ambiguous calls in nanopolish
I am converting scores into frequencies (from nanopolish results) using Nanopore scripts.
The script filters out lots of rows on the premise that reads where LLR < T * N are ambiguous.
where:
LLR = log-likelihood ratio (log-likelihood methylated / log-likelihood unmethylated)
T   = arbitrary threshold. Above that value the C is methylated.
N   = number of Cs (CpG) in the fragment analyzed.

What? Why? Where?

What's the question here? If there is not enough support for a position, it gets filtered out...

I'm assuming that you are using the workflow defined here: https://nanopolish.readthedocs.io/en/latest/quickstart_call_methylation.html . Can you give more information about where you are in the process of the "Nanopore scripts?" I am guessing that you are at calculate_methylation_frequency.py based on available information, and in that case it is as simple as Wouter de Coster says.

Assuming that it is the postprocessing script calculate_methylation_frequency.py that you are using, user Wouter de Coster's comment is correct that the methylation event is filtered out at the frequency estimation stage if it has a low LLR. It will have  a low LLR if it is considered not confident enough to report. This means not only LLR > T * N, but also the statistical threshold LLR > 0 (see below).
Specifically, here is the code in the script:

def update_call_stats(key, num_called_cpg_sites, is_methylated, sequence):
    if key not in sites:
        sites[key] = SiteStats(num_called_cpg_sites, sequence)

    sites[key].num_reads += 1
    sites[key].called_sites += num_called_cpg_sites
    if is_methylated > 0:
        sites[key].called_sites_methylated += num_called_cpg_sites

That is the function that decides whether to add a site to output, which depends on is_methylated, which is set here:
    # [skip several lines of code]
        num_sites = int(record['num_motifs'])
        llr = float(record['log_lik_ratio'])

        # Skip ambiguous call
        if abs(llr) < args.call_threshold * num_sites:
            # [MP NOTE: SKIP if LLR < T * N]
            continue
        sequence = record['sequence']

        # [MP NOTE: SKIP IF LLR <= 0]
        is_methylated = llr > 0

        # [skip several lines of code]

        update_call_stats(key, num_sites, is_methylated, sequence)

It's a little messy to read because the sites dictionary is a global outside of update_call_stats().
So overall, it seems like if you want sites included that are below this threshold, you will have to either:

use an earlier stage of processing (output of nanopolish call-methylation I believe according to the docs) or
alter the script to change the heuristic thresholds.

Thanks! It wasn't easy to follow the script. I understood that the whole fragment is discarded if it can be defined which CpG is methylated (if at all). @Wouter de Coster, is that correct?


text:
Get a comprehensive list of virulence factors for certain or all pathogens?
Is there a database with all known proteins/genes of specific pathogens such as Staph aureus and annotated whether it is a known virulence or resistance factor that idealy can be downloaded in table format?

Please consider accepting MP's answer. For my own stuff I found it really useful, its a good a question and really good answer.

I suggest that you are looking for a resource like VFDB (virulence factor database). That table does annotate by bacterial species, so you can filter as you like.
You can download a description table or FASTA sequences from their downloads page (see screenshot). Last update was 4 days ago, so they should be quite current.

For resistance genes, you could use something like CARD (Comprehensive Antibiotic Resistance Database), which also has a downloads page where you can download a similar TSV table or FASTA sequences.
Both resources include web tools for predicting the presence of the elements of interest in sequences.
If you are looking for a state of the art prediction tool for such elements in genomes, you might look into seqscreen from Todd Treangen's lab.
For a general purpose microbiological annotation and analysis resource, you might also look into BV-BRG, which succeeds PATRIC. I have only used PATRIC but BV-BRG seems very thorough.
I think that you are unlikely to find a resource that has non-virulent/non-resistance genes annotated.

Wow what a cool answer. Got to be the accepted answer, but lets wait. There are a number of microbiological bioinformaticians here.

Thanks for the input. I saw the VFDB database, but I had concers how up to date it is. I am going to check out the other resources you mentioned.

@Soren added a note- VFDB was last updated 4 days ago. See screenshot.


text:
dog coordinates (canFam3) to human coordinates (hg19)
I've converted dog coordinates to human using UCSC LiftOver. These are 200bp intergenic regions that are differentially methylated from normal dogs to cancer dogs. I've converted these to human coordinates and found that a lot of them overlap with differentially methylated regions we found in the human model.

Is this okay to do?
What stringencies should I check or modify, the min ratio of bases needed to overlap?
Lastly, is LiftOver taking the 200bp dog sequence and saying that at least 10% of the sequence has to align with a region within the human genome?

I could not find an answer online.

Based on this old answer, this is not recommended for liftOver. For between-species liftovers, it is claimed that you are likely to want a different tool such as pslMap. For a lot of rather detailed discussion of how liftOver operates in this context and how to parameterize properly, see here.
See also this old Biostars question addressing the issue of mapping not between versions of the same assemblies but different genomes entirely. This thread also notes that using liftOver can be dangerous. That isn't to say that you can't learn something interesting, just that you might make mistakes.
pslMap apparently drops non-syntenic regions, which may lose a lot of data, but it is supposedly making rather fewer assumptions about assembly relationships.
orthoMap is apparently a similar tool to pslMap, but I can't find a way to run it except for a very old and restricted web interface. There are other tools with very similar names ("orthomap") but somewhat different uses, which may increase confusion.
A counterpoint to all the foregoing is that this rather detailed sum-up of liftover tools does not seem to be very concerned about the issues of using tools like liftOver between species. It specifically uses dog to human liftover as the example, just as you want to. It's notable that the assumption is made in that page that a high-quality chain file exists, which does seem to be true for you.
To go through your questions:

depends, see above.
For suggested parameters for running blat/chain generation/liftOver for this application, see here.
This is derived in part from the [chain file] I believe rather than liftOver parameters (http://genomewiki.ucsc.edu/index.php/LiftOver_Howto). However overall this is set by the parameters you used in generating a chain file, which probably boils down to parameters of the workflow here. liftOver itself is just relying on the chain alignments for each region with the best score.


text:
Seeking Guidance on Designing DNA Sequences with Specific Criteria
Dear Bioinformatics Community,
I hope you're all doing well. I am currently working on a project that requires designing some 10 nucleotide single-stranded DNAs, (preferably), with specific criteria.
The criteria for the DNA sequences are as follows:

An equal mix of nucleotides (A, C, G, and T)
Thermodynamically stable sequences
Sequences that do not self-hybridize or hybridize with each other

I have been exploring the NUPACK software suite to achieve these goals, but I am finding it challenging to use the tools effectively for this specific task. I understand that NUPACK is primarily designed for generating sequences based on target secondary structures, which may not directly align with my requirements.
I would greatly appreciate any advice, suggestions, or alternative methods you can provide to help me design these DNA sequences. Are there any other tools or resources that might be better suited for this purpose?
Thank you in advance for your insights and expertise. Your guidance will be invaluable as I work to complete this project.
Best regards,

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

Welcome @Hashem and thank you for your question. It would useful to understand the downstream application, i.e. why do you need the synthetic oligonucleotides and some understanding of the organism, possibly in general terms, e.g. just kingdom level and possible diversity. This appears to be primer design, for which there are a lot packages, within that you appear to be wanting nucleotide degeneracy.  The question is not hard to answer - but is not very clear.

I edited your question to make it clearer and to add relevant tags, please update if you do not consider my edits to be accurate.

It is not difficult to simulate every 10mer DNA oligo, here is some Python code that will do that:
(there is probably a more elegant way to do this)
from itertools import product

nts = ["A", "T", "C", "G"]
# create a tuple with 10 copies of nts for the cartesian product
ten_tup = (nts,)*10

decamers = []
# create and iterate over a stringified cartesian product,
# which consists of all 10mers.
for a,b,c,d,e,f,g,h,i,j in product(*ten_tup):
  decamers.append(
    a+b+c+d+e+f+g+h+i+j
  )

At this point, you have some choices to make:

how many 10mers do you want?
what exactly is your criterion for GC content?
what exactly is your criterion for thermodynamic interactions?

I'd suggest either primer3 or viennarna for evaluating thermodynamics. You can run each from the command line or from a web server, but they do have Python APIs, which I use below. ViennaRNA is definitely more trustworthy for sequences as short as 10, so I'll use it here. It will be a lot slower though as it's not cheating by using the nearest-neighbor approximations.
I have made up some numbers for the thresholds below, replace with your own.
from random import shuffle
GC_MIN = 4
GC_MAX = 6
DDG = -5
NUM_WANTED = 10

def test_gc(seq, gc_min=GC_MIN, gc_max=GC_MAX):
  gc_count = seq.count("G") + seq.count("C")
  if gc_count >= gc_min and gc_count <= gc_max:
    return True
  else:
    return False

def test_thermo(seq, other_seqs, ddg=DDG):
  intra = RNA.fold_compound(seq).mfe()[1]
  print(intra)
  if other_seqs:
    # find the most favorable interaction with any other sequence
    inter = min(
      [RNA.fold_compound(seq+"&"+other_seq).mfe_dimer()[1] for other_seq in other_seqs]
    )
  else:
    inter = 0
  homo = RNA.fold_compound(seq+"&"+seq).mfe_dimer()[1]
  if intra > ddg and inter > ddg and homo > ddg:
    return True
  else:
    print("thermo reject!!")
    return False

chosen_seqs = []
# shuffle the decamers to ensure that you get dissimilar sequences
shuffle(decamers)
i = 0
for seq in decamers:
  i+=1
  print(i)
  if test_gc(seq) and test_thermo(seq, chosen_seqs):
    chosen_seqs.append(seq)
  if len(chosen_seqs) >= NUM_WANTED:
    break

# chosen_seqs has your decamers that pass the threshold

In principle, you could even combine these 2 steps but I think it's cleaner to keep them apart unless storing 4^10 ~= 1 million decamers in memory is an issue, which it shouldn't be with modern computers.
Caveat emptor, this code may contain bugs and I haven't tested it other than making sure it runs.


text:
How to compare CDS, CNE, miRNA, UCNE across new related genomes?
I have 5 sets of raw Illumina whole genome shotgun data from 5 different species, and 3 reference genomes that are assembled/annotated, I want to produce a similar table to Table S2 of this paper: https://www.sciencedirect.com/science/article/pii/S0960982221011908#app2.
This involves comparing reads from each of my 5 new genomes to each of the 3 reference genomes.
How can I go about aligning the sets of seq data to the references so I can capture the CDS, conserved noncoding elements (CNEs), miRNA, and ultraconserved noncoding elements (UCNE)?
I only have CDS annotations for the 3 reference genomes, so I will also need to annotate the miRNA, CNE, and UCNE myself. Suggestions for software to do this would be helpful.

How do the reference genomes map onto the raw data you have? E.g. you have a big species complex or genus with several species, and there are 8 species, of which 3 have reference genomes? Or are some of your sequence data different isolates of the same species that you have references for? THis changes how you want to analyze the data significantly. Furthermore: what is the raw data? Is it Illumina genomic, PacBio genomic, Illumina RNA-seq, PB RNA-seq....?

Ok looking at the table: you have 5 species, each of which you want to compare to each of 3 reference genomes- is that accurate? Also: do you have annotation tracks for the stuff you're interested in on each of your reference genomes?

@MaximilianPress The raw data is Illumina. Yes i do have 5 species, that I want to compare to the 3 reference genomes, im working with lizard species, and the three reference are gekko japonicus, leopard gecko, and crested gecko. The CDS is annotated for each reference, but I believe there isn't miRNA info avaliable.


text:
Error in Load Data in GSEA 4.3.3
I've made my files for a GSEA ANALYSIS, but the software gives me this error:
---- Full Error Message ----
There were errors: ERROR(S) #:1
Parsing trouble
java.lang.NumberFormatException: ...

---- Stack Trace ----
# of exceptions: 1
------For input string: "description"------
java.lang.NumberFormatException: For input string: "description"
    at java.base/jdk.internal.math.FloatingDecimal.readJavaFormatString(Unknown Source)
    at java.base/jdk.internal.math.FloatingDecimal.parseFloat(Unknown Source)
    at java.base/java.lang.Float.parseFloat(Unknown Source)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseStringToFloat(AbstractParser.java:250)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseFieldsIntoFloatArray(AbstractParser.java:361)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.AbstractParser.parseTextMatrixToDataset(AbstractParser.java:277)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.TxtDatasetParser.parse(TxtDatasetParser.java:108)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.readDatasetTXT(ParserFactory.java:200)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:758)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:735)
    at org.gsea_msigdb.gsea/edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:53)
    at java.desktop/javax.swing.SwingWorker$1.call(Unknown Source)
at java.base/java.util.concurrent.FutureTask.run(Unknown Source)
at java.desktop/javax.swing.SwingWorker.run(Unknown Source)
at java.base/java.util.concurrent.ThreadPoolExecutor.runWorker(Unknown Source)
at java.base/java.util.concurrent.ThreadPoolExecutor$Worker.run(Unknown Source)
    at java.base/java.lang.Thread.run(Unknown Source)

These are the first lines of my file, made up by 18839 rows of genes:
#1.2
18839   9
NAME    description ASPC1   CFPAC1  HS766T  BXPC3   CAPAN1  HPAFII  PANC0813    PaTu8902    SU8686
A1BG    NA  0.13750352  0.15055968  0.16349873  0   0.37851162  0.01435529  0   0.01435529  2.10433666
A1CF    NA  3.45285896  0.08406426  0.05658353  0   0.32192809  2.1210154   0.01435529  0   0.04264434
A2M NA  0.07038933  0.20163386  0.48542683  0.21412481  0.05658353  0.04264434  0.20163386  0.05658353  0.15055968
A2ML1   NA  0.01435529  0.01435529  0   1.88752527  0.21412481  0.02856915  2.70708299  2.63459327  0.71369581
A3GALT2 NA  0.02856915  0   0   0   0   0.02856915  0.04264434  0.02856915  0
A4GALT  NA  3.46466827  3.36176836  3.09423607  4.22573765  2.50589093  3.61941301  3.21723072  0.01435529  3.34198575
A4GNT   NA  0   0   0   0   0   0   0   0.01435529  0.26303441
```

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

A similar error has been described in the past and attributed to the presence of missing observations in the data.
Have you checked your data and tried to look for e.g. NA values? It seems that GSEA does not work with missing data.