text:
BLASTn query coverage discrepancy
Using the BLASTn webtool with,
Query Seq: CAATTGCGAAGCCACATTATT
Subject Seq: CP031332.1
BLASTn returns a single hit to CP031332.1 with query coverage at 100%.
However if one looks at the alignment, the query coverage is clearly 17/21 ≈ 81%. Could someone explain to me the reason for this discrepancy?
Edit,
Web version: https://blast.ncbi.nlm.nih.gov/Blast.cgi?PROGRAM=blastn&PAGE_TYPE=BlastSearch&LINK_LOC=blasthome
Output: https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Get&RID=RCK3F55Y114

To expand on the excellent answer by haci:
BLAST reports (a) the summary of alignments and (b) individual alignments, see below. Note that your link to the results expired (as NCBI BLAST results links tend to do). It is considered best practice to post the text version of the relevant portion of the results in the question itself.
The summary provides an overview of all the query sequence individual hits in the subject sequence. Because the individual hits (alignments, or ranges) together cover 100% of the query and each of these has the max identity = 100%, that is what it reports in the summary. It is not very intuitive for many users.
You should always look at the individual hits in detail. The above 100% coverage in the summary consists of ranges 1 and 14, shown below. The top hit (range 1) covers only 17 out of 21 nt with 100% identity, and this is the most relevant information from BLAST. The other hit (range 14) covers 10 out of 21 nt with 100% identity. This is the cause of the discrepancy.
These 2 alignments are not close together in the subject sequence. Depending on your goal, choose one or more of the individual alignments and consider alignment length and % identity of those hits individually, rather than summarized together as in the BLAST report summary. For most applications, the top hit (range 1) would be the most relevant alignment.
Program: BLASTN
Subject:Neisseria meningitidis strain M22814 chromosome, complete genome ID: CP031332.1(nucleic acid) Length: 2190201
Query #1: seq1 Query ID: lcl|Query_56011 Length: 21

Sequences producing significant alignments:
                                                                  Max    Total Query   E   Per.
Description                                                       Score  Score cover Value Ident  Accession
Neisseria meningitidis strain M22814 chromosome, complete genome  34.2   57729 100%  0.001 100.00 CP031332.1

Alignments:
>Neisseria meningitidis strain M22814 chromosome, complete genome
Sequence ID: CP031332.1 Length: 2190201
Range 1: 2074883 to 2074899

Score:34.2 bits(17), Expect:0.001,
Identities:17/17(100%),  Gaps:0/17(0%), Strand: Plus/Plus

Query  1        CAATTGCGAAGCCACAT  17
                |||||||||||||||||
Sbjct  2074883  CAATTGCGAAGCCACAT  2074899

...

Range 14: 354292 to 354301

Score:20.3 bits(10), Expect:15,
Identities:10/10(100%),  Gaps:0/10(0%), Strand: Plus/Minus

Query  12      CCACATTATT  21
               ||||||||||
Sbjct  354301  CCACATTATT  354292

Your blast result is not composed of a single alignment of all of the length of your query vs the subject but from a series of alignments: in each a subsequence of your query is used and these alignments are sorted by "E value", "score", "percent identity", ...
Now in your case, you are talking about the first hit, in which a subsequence of your query (positions 1-17) match perfectly to a portion of the subject. Query coverage is not calculated from this first hit though, it takes all of the alignments into account and again in your case all 21 bases were used in one or more alignments. I believe it would be fair to filter out all but the first hit given the E-values and re-calculate the query coverage as 17/21.


text:
Does rRNA depletion protocol give higher number of mapped reads in Intronic regions?
Recently, I have downloaded a publicly available dataset, which are 350 tumor samples. I see the following information from the published paper.

They used Ribo Zero Gold and rRNA was depleted. Strand specific data. After aligning the data I did some alignment quality check with Qualimap RNA-Seq QC tool. I visualised the bam files in IGV. Alignment is good. For all samples 90% alignment rate was seen. I observed that in all samples Higher percentage of mapped reads were originating in Intronic regions. Followed by Exonic and intergenic regions.
I have seen a post here Reads mapped to exonic, intronic and intergenic regions where they say high intronic reads could be because of contamination. I googled about higher reads in intronic regions and found some papers Evaluation of two main RNA-seq approaches for gene quantification in clinical RNA sequencing: polyA+ selection versus rRNA depletion and some other links RIBO-DEPLETION IN RNA-SEQ – WHICH RIBOSOMAL RNA DEPLETION METHOD WORKS BEST? in which they said Greater intronic reads were with rRNA depletion protocol.
And even in this RNA-SEQ tutorial, it is mentioned that - A higher intronic mapping rate is expected for rRNA removal compared to polyA selection.
So, my question:
I am working with lncRNAs. So, I'm using the samples prepared with rRNA depletion protocol. Is this higher intronic rate is common in rRNA depleted dataset or do I have to check anything else to proceed further with these samples?

Hi @maven .. cool username. I'm not really sure this is a bioinformatics question.

It's not so much that you have "intronic contamination" or "genomic contamination", rather you're not selecting explicitly for full-length mature transcripts with rRNA depletion. That is the most common cause for higher intronic read rates. There's nothing you can do about this post-hoc, just continue along.
BTW, many lncRNA's are polyadenylated, so you'll keep them with poly-A selection.

Thanks for the answer Devon. I saw this - Intronic and intergenic mapped reads could be long non coding RNAs, most of them are unannotated [https://www.researchgate.net/post/What_are_the_causes_of_reads_mapped_to_intergenic_region_in_RNA-seq]. My approach is to look for novel lncRNAs.

Yes, it's a non-trivial problem if you can't filter out immature sequences before library prep.

Just to make sure, I also checked the rRNA contamination. And I observed there is very less percentage. 12408836 + 0 in total (QC-passed reads + QC-failed reads) 0 + 0 secondary 0 + 0 supplementary 0 + 0 duplicates 81734 + 0 mapped (0.66% : N/A) 0 + 0 paired in sequencing 0 + 0 read1 0 + 0 read2 0 + 0 properly paired (N/A : N/A) 0 + 0 with itself and mate mapped 0 + 0 singletons (N/A : N/A) 0 + 0 with mate mapped to a different chr 0 + 0 with mate mapped to a different chr (mapQ>=5)

The read counts I extracted using featureCounts are low compared to total alignments number. Is this because the read counts are from the Exonic regions?

Yes, featureCounts ignore intronic reads by default.

Are there any tools to check DNA contamination?

Probably, you’d want to see how much intergenic expression you have, though you’d need some baseline for comparison.

For all the 350 samples I downloaded...I see 15-30% exonic reads, 70-80% intronic reads an 5-10% interegic reads. so, what % of intergenic tell us there is DNA contamination?

You’d need to mine the literature, but that’s not unreasonable levels given that most of it is probably rRNA and tRNA.

Hi Devon....How to filter our rRNA from the data....And you said more intronic reads could be from immature transcripts....Is it possible to use such data for novel lncRNA detection?

It would be very difficult to use this data for lncRNA detection if you can't be relatively certain that it's high quality. You're not going to have an easy go of it. Regarding rRNA exclusion, blacklist very high expression regions that end up having tRNA or rRNA homology.

I have small basic question...How do I get the sequencing depth of the dataset I downloaded.

Count the number of reads and look at their length distribution. FastQC can aid in this.

I can't possibly see how intron contamination is linked to removal of rRNA depletion.
The only reason it would appear to have increased the number of contaminants is because post-rRNA removal the proporation of intron contaminants has increased against the total remaining RNA content. However, the actual total number of contaminants remains exactly the same pre- and post-rRNA removal. By the same token the RNA of interest lncRNAs will also have proportionally increased, so you get a better depth of its predominance and diversity.
Thats just life and perhaps just filter this bioinformatically.

Hi @Michael thanks for the answer. You mean I should remove these rRNA contamination and then proceed? How do I do that? Any tips please

source: Stack Exchange

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

sort: votes


text:
How to get rows with similar values in two different columns using command line?
I have following example data from a large file.
cont100 1128 1125
cont1005 3642 3642
cont1006 103 19
cont1037 3146 3146
cont104 895 890
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577
cont1078 43 42
cont1081 727 2

I need to get only the rows with similar values in column 2 and column 3. The output should be like this.
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Also I need to get the inverse of that. That means I need to get the rows with values that are not similar in the columns. The output should be like this.
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

Do you mean identical values in columns 2 and 3?

I think it is easy to do with awk, e.g. see here.
For your example, with your data in file:
%  awk '$2 != $3 {print $0}' file
# outputs
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

%  awk '$2 == $3 {print $0}' file
# outputs:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Note that where you say "similar" I am interpreting that to mean "exactly the same", because that is the case in your example.

It worked and I got the result that I wanted. Thank you very much.

The Perl solution is similar to the awk solution from MaximilianPress:
Print lines with identical values in columns 2 and 3:
perl -lane 'print if $F[1] == $F[2];' in.txt

# Prints:
cont1005 3642 3642
cont1037 3146 3146
cont1056 934 934
cont1059 1750 1750
cont1072 2577 2577

Print lines with different values in columns 2 and 3:
perl -lane 'print if $F[1] != $F[2];' in.txt

# Prints:
cont100 1128 1125
cont1006 103 19
cont104 895 890
cont1078 43 42
cont1081 727 2

The Perl one-liner uses these command line flags:
-e : Tells Perl to look for code in-line, instead of in a file.
-n : Loop over the input one line at a time, assigning it to $_ by default.
-l : Strip the input line separator ("\n" on *NIX by default) before executing the code in-line, and append it when printing.
-a : Split $_ into array @F on whitespace or on the regex specified in -F option. Arrays in Perl are 0-indexed.
SEE ALSO:
perldoc perlrun: how to execute the Perl interpreter: command line switches

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


text:
Easy way to create gene plot in gviz?
I would like to create a simple plot of a given human gene containing the main features, e.g. exons, introns, ORFs, transcripts etc. I understand that gviz  is a popular choice for this kind of thing, but most tutorials are massive overkill for what I want to do (I don't need fancy features such as importing BAM files and plotting coverage).
I was wondering whether someone can point me towards some sample code for creating a basic gene plot. It probably would be good to start by creating a dataframe containing all the features of my gene - is there a handy package that can do that?

If Python is an option, you could check BioPython and GenomDiagram. This should get the job done.

I assume you are trying to make something like a "gene arrow map" and [gggenes][1] can help you (with some little tweaks) in case you are familiar with ggplot2.


text:
Error while using E-utilities on slurm HPC
I keep getting the error while using E-utilities on Slurm HPC, could you suggest how I can resolve it?
SCRIPT:
#!/bin/bash
#SBATCH --job-name=tax_job_test      # Job name
#SBATCH --mail-type=END,FAIL         # Mail events (NONE, BEGIN, END, FAIL, ALL)
#SBATCH --mail-user=na@gmail.com    # Where to send mail
#SBATCH --ntasks=8                  # Number of MPI ranks
#SBATCH --nodes=1                    # Number of nodes
#SBATCH --time=23:05:00              # Time limit hrs:min:sec
#SBATCH --partition=k2-medpri

module load e-utilities/03.02.19
module load apps/perl/5.20.2/gcc-4.8.5

for i in `cat taxonomy.id.txt` ; do efetch -db nuccore -format gbc -id ${i} |xtract  -pattern INSDSeq -element INSDSeq_accession-version INSDSeq_taxonomy INSDSeq_organism >>taxo.2020.txt; done

input taxonomy.id.txt is like this:
AKN43077.1
AKN41814.1
PNE10086.1
AKN42672.1
AKN41654.2
AKN43077.1
PNE10086.1
AKN43189.1
TMK08501.1
TMK08501.1
AKN42976.1
AKN43189.1
AIZ47938.1

Error:
The script keep giving below error:
501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed)
No do_post output returned from 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network'

Result of do_post http request is

$VAR1 = bless( {
                 '_rc' => 501,
                 '_request' => bless( {
                                        '_method' => 'POST',
                                        '_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
                                        '_content' => 'db=nuccore&id=AKN43077.1
&rettype=gbc&retmode=xml&edirect_os=linux&edirect=11.7&tool=edirect&email=3052771@node126.pri.kelvin2.alces.network',
                                        '_headers' => bless( {
                                                               'content-type' => 'application/x-www-form-urlencoded',
                                                               'user-agent' => 'libwww-perl/6.39'
                                                             }, 'HTTP::Headers' )
                                      }, 'HTTP::Request' ),
                 '_content' => 'LWP will support https URLs if the LWP::Protocol::https module
is installed.
',
                 '_headers' => bless( {
                                        'client-warning' => 'Internal response',
                                        '::std_case' => {
                                                          'client-warning' => 'Client-Warning',
                                                          'client-date' => 'Client-Date'
                                                        },
                                        'client-date' => 'Mon, 05 Oct 2020 14:27:37 GMT',
                                        'content-type' => 'text/plain'
                                      }, 'HTTP::Headers' ),

Many thanks

501 Protocol scheme 'https' is not supported (LWP::Protocol::https not installed) indicates that a specific perl module is missing.
If indeed your edirect installation is from March 2019 (as I interpret from e-utilities/03.02.19), you should update your edirect installation to the latest version. If you update, I expect the necessary perl module will be installed and the error will go away.

Ah! of course ... consider using CPAN to load the module

I just downloaded the script install-edirect.sh from this link ((https://www.ncbi.nlm.nih.gov/books/NBK179288/), then execute it by running:
source ./install-edirect.sh
and it works well!
Many thanks


text:
How to identify unknow bacterium species from whole genome genetic sequence
I am a Biochemist that is unfamiliar with bioinformatic tools and new to academia as a whole.
I am currently using ILLUMINA PE data, which I trimmed (Trimmomatic), corrected (Rcorrector) and assembled (SPAdes). I am now interested in using the genetic sequences from my contigs to identify the original organism (I know it is a bacterium).
I have tried using BLASTn (with default parameters) on one of my contigs. I let the query run for 30 minutes before stopping, nothing came of it.
What would be the best way to go about doing this?

posted an answer here for the cross-post. https://biology.stackexchange.com/questions/96183/how-to-identify-an-unknown-species-from-a-fasta-file-containing-its-genome-sequ/96198#96198

I think the best option is using this online tool (https://blast.ncbi.nlm.nih.gov/Blast.cgi) If you don't have too many contigs you can just do that manually. It usually takes less than one minute even with large contigs. But if it takes too long... just use a shorter sequence and it will be as effective to find the organism.

You can use KmerFinder and Kraken2 for taxonomy assignment


text:
Seurat Violin Plot: Why do dots align in one row?
I am analyzing a rather small scRNAseq dataset (39797 features across 3522 samples) and I have observed that dots in the violin plot are plotted in one line next to each other rather than in a cloud of dots (see plot below):

I am just wonndering what the reason might be for this?
I have followed the Integration and Label Transfer tutorial using SCTransform (https://satijalab.org/seurat/v3.1/integration.html). After the integration process I switched back to the SCT assay, which should contain the normalized counts. For plotting I used the following command:
VlnPlot(..., assay = "SCT", slot = "data"), again specifying to use the normalized counts.
Can someone explain me why this is the case?
Thank you already in advance!

That is because SCTransform (mostly) outputs discrete values, you can check the @slot of your Seurat object.

Thank you for your answer!Do you know why this is mostly for lower expression values?

It's log-scaled.


text:
Finding motifs: fasta file with 10,000 sequences
I am new to Python. I am trying to parse a fasta file containing 10,000 sequences to look for motifs (microsatellites in particular). I tried using Seq Utils to parse my sequences for a particular motif, that is, (TG)20. I suspect because I used the string function (str), it treated the entire fasta file as a one single fasta sequence. I then tried using BruteForce on the motif (TG)17 and found that an error appeared which said that "...the object of type 'FastaIterator' has no len()...". The script however worked when only a single sequence was in the fasta file.
In essence, I need some help figuring out how to:

Prepare an instances object(is this correct?) with microsatellite
motifs (for starters, I intend to parse the sequences for dinucleotide repeats (TG)10,(AC)10,(AT)10
and (AG)10;
Parse my 10,000 long sequence file using this object(?); and
Print an output specifying which of my sequences contain microsatellite motifs.

Appreciate any help I can get on this.
#Seq Utils to find Motifs
>>> pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
>>> sequence = Seq(a)
>>> results = SeqUtils.nt_search(str(a),pattern)
>>> print(results)
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 38, 40, 42, 44, 46, 48, 50, 52, 54,
 56, 58, 68682, 68684, 68686, 370986, 370988, 370990, 370992, 370994, 370996, 43
8620, 438622, 438624, 438626, 438628, 438630, 438632, 438634, 438636, 438638, 43
8640, 556703, 556705, 556707, 784973, 784.... # it's too long so I am just pasting the first few lines

#BrutForce to find Motifs in a fasta file containing 10,000 sequences
>>> def BruteForce(s, t):
...     occurrences = []
...     for i in range(len(s)-len(t)+1):
...         match = True
...         for j in range(len(t)):
...             if s[i+j] != t[j]:
...                 match = False
...                 break
...         if match:
...             occurrences.append(i)
...
...     print(occurrences)
...
...

>>> t = "TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG"
>>> len(t)
34
>>> s = SeqIO.parse("group23.fas", "fasta")

>>> BruteForce(s,t)
Traceback (most recent call last):
  File "<console>", line 1, in <module>
  File "<console>", line 3, in BruteForce
TypeError: object of type 'FastaIterator' has no len()

#BrutForce to find Motifs in a fasta file containing a single sequence
>>> s = SeqIO.read("group23v.fas", "fasta")
>>> BruteForce(s,t)
[0, 2, 4, 6, 8, 10, 12, 14, 16, 18, 20, 22, 24, 26, 28, 30, 32, 34, 36, 38, 40, 42, 44, 46, 48, 50, 52, 54, 56, 58, 60, 62, 64, 66, 68, 70, 79, 98, 117, 175, 193, 198, 217, 221, 277, 279, 282, 289, 294, 304, 311, 335, 381, 437, 464, 468, 510, 514, 522, 593, 602, 622, 634, 647, 658, 678, 704, 727, 750, 773, 804, 889, 895, 908, 911, 922, 937, 969, 975, 978, 1002]

Bio.SeqIO.parse() returns a SeqRecord iterator. To access the sequence of this sequence record, you need to use the .seq attribute so you should update your s with s.seq.

To expand on haci's answer (and since you mentioned being new to Python) you can loop over the iterator object and use the seq attribute like this:
pattern = Seq("TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG")
for record in SeqIO.parse("group23.fas", "fasta"):
    results = SeqUtils.nt_search(str(record.seq), pattern)
    print(results)

(More info on iterators here.)  That will work the same whatever the number of sequences in group23.fas.
I'd be cautious about looking for microsatellites with fixed strings.  If the sequences you're looking for are short and have imperfections among the repeats, you might miss them, depending on how many perfect repeats your search pattern has.  For example with these stub example sequences:
>seq1
TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq2
ACGATGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG
>seq3
TGTGTGTGTGTGTGTGTGTGTCTGTGTGTGTGTGTGTGTG

The nt_search output gives you:
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 0]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG', 4]
['TGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTGTG']

That is, one match at position 0 of the first sequence, one at position 4 of the second, and no matches for the third (there's a C I put in the middle).  This might not be relevant for you if you're looking microsatellites much longer than your search patterns, but for shorter ones it might be a problem.

Hi Jesse and @Haci. Thanks so much! I followed both your advice and it worked. I tried something less common, that is, trinucleotide repeats (AGC)10 i.e. pattern=Seq('AGCAGCAGCAGCAGCAGCAGCAGCAGCAGC')
I used the following script

for record in SeqIO.parse("group23.fas", "fasta"):
...     results = SeqUtils.nt_search(str(record.seq), pattern)
...     print(results)
...     print(record.id)

The nt_search output gave the following! (only showing two examples):
NODE_226753_..cov_2.031008 ['AGC...']
NODE_226754_..._cov_2.021041
['AGC...', 265, 268, 271, 274]

Here's a solution that I came up with after help from both haci and Jesse. It involves using the following modules: os, glob, SeqIO, Seq, SeqUtils.
import os, glob

os.chdir("C:\python38\Lib\site-packages\Bio\SeqIO")

from Bio import SeqIO

from Bio.Seq import Seq

from Bio import SeqUtils

patttern = Seq("ACACACACACACACACACAC")

folder_path = 'C:\python38\Lib\site-packages\Bio\SeqIO' #specified the folder where my sequences were stored

fasta_paths = glob.glob(os.path.join(folder_path, '*.fasta')) #glob module finds all the pathnames matching a specified pattern i.e. '*.fasta'

with open("motif4.txt", "w") as f:
... for fasta_path in fasta_paths:
...     print(fasta_paths)
...     for seq_record in SeqIO.parse(fasta_path, "fasta"):
...         results = SeqUtils.nt_search(str(seq_record.seq),pattern)
...         f.write(str(results)) #replace with seq_record.seq to save the list of node names of all sequences that were parsed
...         f.write("\n") #ensures results of the parse through all the sequences in the fasta_path is written to a new line
... f.close()

The nt_search output gave the following output (similar to what Jesse showed):
1 ['ACACACACACACACACACAC']
2 ['ACACACACACACACACACAC']
3 ['ACACACACACACACACACAC']
4 ['ACACACACACACACACACAC']
5 ['ACACACACACACACACACAC']
6 ['ACACACACACACACACACAC', 16636, 16638, 16640, 16642, 16644, 16646, 16648, 16650, 16652, 16654, 16656, 16658, 16660, 16662, 16664]
7 ['ACACACACACACACACACAC']
8 ['ACACACACACACACACACAC']
9 ['ACACACACACACACACACAC']
10 ['ACACACACACACACACACAC']
11 ['ACACACACACACACACACAC', 9174]

It would miss sequences that did not match (AC)10 or those with imperfect repeat motifs. However, in cases where the motif was > (AC)10, nt_search listed down all locations where the motif was found sequentially. For example, in line 6, all locations were listed and showed the motif was in fact 36 nucleotides long in total (pos #16,636 to #16,673).
You can then import the txt. files into a databasing app and analyse the results.


text:
How do I obtain a vcf read from a pyvcf.Reader?
I wondered how would I be able to obtain variant based on chrom and location range.
Optimally, this function would be like:
import vcf
vcf_reader = vcf.Reader(filename=some_file_name)
vcf_read = vcf.variants_between(CHROM, START, END)

Does anyone how to do so?

I think vcf_reader.fetch is what you want.  You have to index the file first though.
If starting from regular text VCF you can run:
bgzip example.vcf
tabix example.vcf.gz

Then the Python code is almost the same as yours:
import vcf
vcf_reader = vcf.Reader(filename="example.vcf.gz")
vcf_read = vcf_reader.fetch("NC_023663.1", 7803980, 7804000)
for record in vcf_read:
    print(record)

For this example VCF you get:
Record(CHROM=NC_023663.1, POS=7803981, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803982, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803983, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803984, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803985, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803986, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803987, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803988, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803989, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803990, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803991, REF=A, ALT=[G, <*>])
Record(CHROM=NC_023663.1, POS=7803992, REF=G, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803993, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803994, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803995, REF=T, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803996, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803997, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803998, REF=C, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7803999, REF=A, ALT=[<*>])
Record(CHROM=NC_023663.1, POS=7804000, REF=T, ALT=[<*>])

API reference:
https://pyvcf.readthedocs.io/en/latest/API.html#vcf.Reader.fetch
Note the "zero-based, half-open coordinate system" it describes (which is why 7803980 in the second Python argument becomes POS=7803981 from the VCF but then it stops at POS=7804000).


text:
How many entries in PDB have SS bonds?
Is there a way to figure out what percentage of entries in PDB have SS bonds?
In addition, what other types of bonds in PDB may be worth looking into (not H bonds)?

The easiest way is to check annotations in mmCIF or PDB files.
Looking at my page with statistics of mmCIF tags, last updated in March:
_struct_conn_type.id        59137×metalc 42706×covale 36195×disulf 10725×hydrog
_struct_conn.conn_type_id   5350402×hydrog 1733277×metalc 559092×covale 226472×disulf

Which means 226472 disulfide bonds in 36195 PDB entries at the time when the PDB had 161470 entries.

Great, thanks! Exactly what I needed.
What about ionic interaction (salt bridges), do you have the numbers?

No, I'm afraid I can't help with salt bridges.


text:
Filtering pileup from site lists
I want to write a script that filters a pileup file from a site lists file.
As an input I get a reference genome, pileup and site lists files.
Example of an output for this script:
Pileup File :
        NC_000001.10    11456   A   0   *   *
        NC_000001.10    11457   C   0   *   *
        NC_000001.10    11458   T   0   *   *
        NC_000002.11    11460   G   1   ,   E
        NC_000002.11    13914   G   1   ,   E
        NC_000003.11    10503990    C   1   ,   E

Reference File (After using grep by '>'):
        >NC_000001.10 Homo sapiens chromosome 1, GRCh37.p13 Primary Assembly
        >NT_113878.1 Homo sapiens chromosome 1 unlocalized genomic contig, GRCh37.p13 Primary Assembly
        >NT_167207.1 Homo sapiens chromosome 1 unlocalized genomic contig, GRCh37.p13 Primary Assembly
        >NC_000002.11 Homo sapiens chromosome 2, GRCh37.p13 Primary Assembly
        >NC_000003.11 Homo sapiens chromosome 3, GRCh37.p13 Primary Assembly

Lists File:
        chr1    11456   C1orf186    -   intronic    intronic    no  no  N   N   N
        chr1    11457   C1orf186    -   intronic    intronic    no  no  N   N   N
        chr2    13914   intergenic  -   intergenic  intergenic  no  no  N   N   N
        chr7    30504355    NOD1    -   intronic    intronic    no  no  N   N   N
        chr3    10503990    SSX2IP  -   Syn Gln->Gln    no  no  N   N   N
        chr1    13131320    MEF2A   +   intronic    intronic    no  no  N   N   N

Output File:
        NC_000001.10    11456   A   0   *   *
        NC_000001.10    11457   C   0   *   *
        NC_000002.11    13914   G   1   ,   E
        NC_000003.11    10503990    C   1   ,   E

So for every refseq in the pileup file I would convert it to the matching chromose in which is found in the reference file, and then overlap it with my lists file so I would get the lines with the same chromose and position.
The problem is that while there are many refseq coordinates that belong to chr1 for example, alignment to each of them will give a position, with respect to the specific coordinate. But on the list there is a full, single, chr1.
I need to be able to compare between:

a pileup file that was created based on a refseq reference genome
a file with editing sites that was downloaded from RADAR (we may give the file)
Maybe there is a more appropriate reference genome file or that there is a way to convert the refseq positions according to the list

EDIT: If I were to convert my pileup file into a bed file (by repeating the position twice) could I use somehow bedtools to solve my problem?


text:
data visualization RNAseq : scaling data for PCA and cluster dendogram
I have count data from a RNAseq experiment (2 samples are from normal cells and 3 samples are cells with a disease), and the data is already standardized by trimmed mean of M values (TMM).
I want to do some plots: biplot of Principal Component Analysis (PCA) and a cluster dendrogram to see if the samples normal/disease are well separated (there is a clear difference between them).
Since the data is already standardized (TMM) should I scale and center the data prior to perform PCA and cluster dendrogram??
thank you!

by TMM you meaning running limma or edgeR on it i supposed? Or obtaining the cpms?

Using kallisto I get the TMM, and before looking for differentially expressed genes with edgeR I am doing PCA and cluster dendogram to look into the data.

Scaling (or centering) makes the genes comparable: Putting the expression levels of genes in the same scale (i.e. between 0 and 1) sustains that all of your genes contribute equally to the PCA or distance calculations. On the other hand, without this step, such calculations would be dominated by the highly/lowly expressed genes. I get "better" results with scaling.

Yes, I understand that. But the data has already been standardized using by trimmed mean of the value. And when performing PCA with and without scaling the results are pretty similar. I am wondering if the scaling and centering are not needed because the data has already been standardized.

What TMM does is to "make different samples comparable" by adjusting for library size and more. It does not deal with individual genes. So yes, scaling would help when expression levels of different genes vary a lot (it is often the case). Since you are using edgeR, you can use the function cpm() to calculate log2 of CPM and then scale these before clustering for example.

ok, thank you for your answer now I undestand. But, in terms of using the cpm() function: is it ok to use it over the data that is already TMM? Thanks


text:
GWAS phenotype data format and preprocessing
I have a set of different phenotypes which I want to use for a GWAS analysis (general linear model). I have a couple of questions and uncertainty about the phenotype data input.

I have control and samples and usually we normalize all sample to the control for other calculations.

Example:
Mean control - 60%
Mean sample1 - 20 %
Mean sample1 normalized to control - 20*100/60 = 33,3 %
Is this normalization step correct for a GWAS analysis?

Should the phenotype data scaled into 0...1 range? I have different scales for different phenotypes but the results should be compared and compatible.

Some of my phenotypes are not normality distributed. What is the best way to handle such data? I have read that the data should be normally distributed for GWAS analysis.

Here are the histograms of two phenotypes:

Please could you clarify if by "GWAS analysis" you are referring to general linear modelling (GLM)? This is regression analysis assessed using regression weights and a residual to understand the association between alleles and phenotypes. Is that correct?

Normalisation is standard procedure for GLM

Yes it's general linear modelling (written in the first sentence).

Okay you need to do this via a Q-Q plot, or at least it is one way to do this. The an alternative is to minimise the residual. When I get a chance I'll post an explanation. You don't necessarily need to transform, it depends on the 'diagnostics', e.g. Q-Q.


text:
alpha diversity wilcox.test
I am trying to do Wilcox test to detect the significant difference in alpha diversity but it is showing an error? physeqN2 is a phyloseq object and Season is a metadata column.
results = estimate_richness(physeqN2, measures = 'Chao1')

d = sample_data(physeqN2)

# calculate Wilcox Test
Wet.root = results[d[,'Season'] == 'Wet',]
Dry.root = results[d[,'Season'] == 'Dry',]
pmv <-pairwise.wilcox.test(Wet.root, Dry.root)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error;
Error in wilcox.test.default(xi, xj, paired = paired, ...) :
  'x' must be numeric
Calls: pairwise.wilcox.test ... compare.levels -> wilcox.test -> wilcox.test.def                                                                                ault
Execution halted

I hvae tried to do it in another way but here also getting error?
alpha.diversity <- estimate_richness(physeqN2, measures=c("Chao1"))

## Wilcox on observed richness
data <- cbind(sample_data(physeqN2), alpha.diversity)
pmv <-pairwise.wilcox.test(Chao1 ~ Season, data)

capture.output(pmv, file = "chao1.16s.wet_vs_Dry.root.txt")

Error:

Error in if (length(ans) == 0L || as.character(ans[[1L]])[1L] == "~") { :
  missing value where TRUE/FALSE needed
Calls: pairwise.wilcox.test ... sapply -> lapply -> FUN -> compare.levels -> [ -> [.formula
Execution halted

Many thanks

Is the data truely paired?

?pairwise.wilcox.test, Your syntax is wrong in both examples pairwise.wilcox.test(results,d) seems more likely. But as Michael mentioned you might not want to do a paired test.

I think this is a bug in the input data.
At a wild guess what I think is happening is the input data is not regarded as paired, hence the error paired = paired. The other weird error is x must be numeric. Firstly,

how are you scoring missing data, its not via - by any chance? Missing data must be via nan type syntax.
are you sure there are no heading bugs? Like a header is getting dragged into the numeric data?

Wilcoxan is a non-parametric T test, so the wet root and dry root must match exactly the same sample. Thus it is not a general wet root versus generic dry root, but wet and dry from the same specimen. If this is the data you have the program doesn't think it is matching.
If it is truely paired, again I think the source of error is you have some missing value which is leaving the pairings unequal lengths or else a header bug.
If the data genuinely isn't paired, then I would look at the Mann-Whitney U test, it would be easier in this scenario.
I don't use R syntax (just pandas dataframes) so can't say whether the syntax in the first example is valid, but it is pretty clear regarding your input. In Python you'd generally just give it the entire dataframe and tell it to go figure.
The reality is that loads of programs do Wilcoxan, you don't need to use R if its a one-off test. Obviously if this a part of a pipeline (automation) thats a bit different.

Hi Michael, I have merged two factor (region & season) into single one and then did annova test followed by posthoc and it work well. Should I update the code?

No that's fine thanks. You can answer your own question however, but are not obliged to do so.


text:
Call NCBI Multiple Sequence Alignment Viewer with target URL?
I am storing a sequence aligment file on a webserver and now would like to visualize the results. I used the NCBI Multiple Sequence Alignment Viewer . The issue is that there are several manual steps involved, that could be streamlined.

open the file dialog
search the add file or type in the urls
click on open

etc.
I wonder if there is a way to call the viewer (or any other alignment viewer) with a single HTTP request, that contains the URL of the alignment file?

Such functionality hasn't been implemented in NCBI's MSA viewer. This is of interest to the developers of MSA viewer and would appreciate feedback from you. Please consider contacting the devs by using the "Feedback" link present at the top right corner of the MSA Viewer.


text:
gene-gene correlation from two different Tissues
I am stuck on how to do correlation for two independent data sets with common row and column names.
A and B are datasets that contain as many rows as genes and as many columns as samples.
The rows in A and B represent a common set of genes but measured in two different tissues.
The columns represent measurements in the same 5 samples in both A and B.
I want to do a correlation between the set of genes in A and B. This is to see if the same genes in both tissues are correlated or not.
Since the matrix would be big in my actual data, I only want to retain a correlation coefficient higher than 0.5.
Here I simulate the data set.
set.seed(1)
A <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(A) <- paste0("G_", 1:100)
colnames(A) <- paste0("M_", 1:5)

set.seed(42)
B <- data.frame(rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100),
                rnorm(100))
row.names(B) <- paste0("G_", 1:100)
colnames(B) <- paste0("I_", 1:5)

Thank you!

You can use mapply(). ta and tb being transposed data frames of your A and B data frames respectively:
> mapply(cor, ta, tb)[mapply(cor, ta, tb) > 0.5]

      G_3       G_5       G_9      G_10      G_11      G_15      G_20      G_23      G_25      G_26      G_33      G_40      G_43      G_48      G_57      G_60
0.5346591 0.8066507 0.8379777 0.6752681 0.7221359 0.5285787 0.7333045 0.5627962 0.6533379 0.7256878 0.5996492 0.6486557 0.5108215 0.7386332 0.6596823 0.6919915
     G_63      G_72      G_76      G_80      G_81      G_90      G_97      G_98      G_99
0.5589583 0.8391917 0.7608801 0.8003665 0.6364557 0.5030968 0.7298439 0.5693024 0.5709411

However, it returns NAs in my studio!
tb = t(B); ta = t(A)mapply(cor, ta, tb)[mapply(cor, ta, tb) > 0.5]

It must be something to do with mapply() treating matrices and data frames differently, not sure. t() outputs matrices not data frames so you can do this: tb = as.data.frame(t(B)).


text:
amino acids substitution probability
I am performing amino acid analysis but require assistance in the underlying amino acid substitution model, i.e. the predetermined rate/frequency of amino acid mutations which are used for example in constructing phylogenetics trees using maximum likelihood or Bayesian methods.
There are numerous matrices of probabilities for every combination of amino acids (something like BLOSUM or PAM matrices and I am wondering within the large number of models which the most appropriate. I am aware there is no 'universal' matrix.
I want to use the matrix to check the probability of amino acid substitution in sequence on a pairwise basis.
For example, I'm not sure whether the BLOSUM or PAM matrix will be useful to me.
[EXTRA] If all diagonal values in such a matrix would be equal to 1 that would be helpful.

The LG model is here:
0.425093
0.276818 0.751878
0.395144 0.123954 5.076149
2.489084 0.534551 0.528768 0.062556
0.969894 2.807908 1.695752 0.523386 0.084808
1.038545 0.363970 0.541712 5.243870 0.003499 4.128591
2.066040 0.390192 1.437645 0.844926 0.569265 0.267959 0.348847
0.358858 2.426601 4.509238 0.927114 0.640543 4.813505 0.423881 0.311484
0.149830 0.126991 0.191503 0.010690 0.320627 0.072854 0.044265 0.008705 0.108882
0.395337 0.301848 0.068427 0.015076 0.594007 0.582457 0.069673 0.044261 0.366317 4.145067
0.536518 6.326067 2.145078 0.282959 0.013266 3.234294 1.807177 0.296636 0.697264 0.159069 0.137500
1.124035 0.484133 0.371004 0.025548 0.893680 1.672569 0.173735 0.139538 0.442472 4.273607 6.312358 0.656604
0.253701 0.052722 0.089525 0.017416 1.105251 0.035855 0.018811 0.089586 0.682139 1.112727 2.592692 0.023918 1.798853
1.177651 0.332533 0.161787 0.394456 0.075382 0.624294 0.419409 0.196961 0.508851 0.078281 0.249060 0.390322 0.099849 0.094464
4.727182 0.858151 4.008358 1.240275 2.784478 1.223828 0.611973 1.739990 0.990012 0.064105 0.182287 0.748683 0.346960 0.361819 1.338132
2.139501 0.578987 2.000679 0.425860 1.143480 1.080136 0.604545 0.129836 0.584262 1.033739 0.302936 1.136863 2.020366 0.165001 0.571468 6.472279
0.180717 0.593607 0.045376 0.029890 0.670128 0.236199 0.077852 0.268491 0.597054 0.111660 0.619632 0.049906 0.696175 2.457121 0.095131 0.248862 0.140825
0.218959 0.314440 0.612025 0.135107 1.165532 0.257336 0.120037 0.054679 5.306834 0.232523 0.299648 0.131932 0.481306 7.803902 0.089613 0.400547 0.245841 3.151815
2.547870 0.170887 0.083688 0.037967 1.959291 0.210332 0.245034 0.076701 0.119013 10.649107 1.702745 0.185202 1.898718 0.654683 0.296501 0.098369 2.188158 0.189510 0.249313

Where the amino acid frequencies are:
0.079066 0.055941 0.041977 0.053052 0.012937 0.040767 0.071586 0.057337 0.022355 0.062157 0.099081 0.064600 0.022951 0.042302 0.044040 0.061197 0.053287 0.012066 0.034155 0.069147
Amino acid order (matrix and frequencies)
A R N D C Q E G H I L K M F P S T W Y V
Reference
Le S.Q., Gascuel O. Molecular Biology and Evolution. 2008 25(7):1307-20.


text:
Looking for genes that are expressed the same in all the samples
I have 5 samples (from a RNA seq experiment) from patients but I have no control, so since we can't look for the differentially expressed genes we want to look for the genes that have similar expression in all the sample.
We want to do that to see if some genes that can be involved are among those expressed in all the samples.
I have never done that, I do not know how to do it and I cant find any info about it, so any ideas are welcome.
Thanks!

Can you be a bit more specific? The vast majority of genes will probably be similar, what question do you want to answer?

Genes "that are expressed the same" would be "the genes that are NOT differentially expressed" so if you cannot study the latter you cannot study the former.

Please [edit] your question and give us some details. Do you already have expression data? What data do you have? What kind of format are they in? Why would you look at all genes that are similarly expressed in all samples? That won't help you identify genes related to a phenotype, not without a control.

standard deviation/mean?

I have edit the question, I hope now is more clear. thanks!

Uhm... if your patients are humans I guess you can find controls on the Internet.

If they are not humans (!!) then you don't need differential expression analysis. Just make expression analysis and then make a script to select transcripts that are equally expressed. It's only a few lines in any language. I would use python for that.

Patients are humans but It is not easy to find controls on the Internet since the protocol and the sequencing has to be done exactly the same way,

I cannot think of a way without performing differential expression analysis. Since you don't have controls, what you could do is assigning some of the samples as controls.
Let's say you have samples A, B, C, D, E and F. You can perform differential expression analysis A, B, C (i.e. treatment) versus D, E (control). If you perform this for all possible combinations (that would be 5 choose 2 which is equal to 5 choose 3 if you decide to go for 3 vs 2 splits) and filter for genes which do not show up in any DE list you end up with the genes that are expressed similarly in al of these samples.


text:
Problems replicating Kover paper results
I am trying to replicate the results for the creator of this repo:  https://github.com/aldro61/kover2_paper for genotype to phenotype machine learning interpretable (decision tree) prediction.
Specifically in the data section when I try to create the dataset as explained here:
https://github.com/aldro61/kover2_paper/tree/master/data
I make a directory called genomes_dir in the data directory. Then I do:
python download_genomes.py 'mycobacterium tuberculosis/kanamycin/metadata.tsv' genomes_dir

at the end of the download it tells me that it saved file genome_paths.tsv the data directory, not genomes_dir. So I run command:
kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress

and it just hangs with a prompt:
(kover) nobu@gold3forever:~/Desktop/BioInformatics/ResearchAssistant/kover2_paper/data$ kover dataset create from-tsv --genomic-data genome_paths.tsv ---phenotype-description "Kanamycin resistance" --phenotype-metadata 'mycobacterium tuberculosis/kanamycin/metadata.tsv  --output data_example.kover --progress
> ^C

it does nothing after hours and hours so I control-C'd to end the process.
Does anyone have any idea what is happening?
I also had to install this to get the datasets from PATRIC:  https://github.com/aldro61/patric_tools

The single quote at 'mycobacterium is not matched so the command prompt is waiting indefinitely for the command to be completed. Maybe you want 'mycobacterium tuberculosis/kanamycin/metadata.tsv' (note closing quote)

Thanks for the catch but now I'm getting some kmer_count = get_kmer_count(tsv_path)   File "/home/nobu/anaconda3/envs/kover/lib/python2.7/site-packages/kover/dataset/create.py", line 136, in get_kmer_count     raise Exception() Exception so I'd better check the data I downloaded.


text:
smallest length of reads in 3'quantseq?
I am studying about RNA seq, and especially about 3'Quantseq(tagseq, 3primeseq). I wonder if there is a cutoff for the reads length. By this I mean that given that 3'quantseq targets the end of the genes and 3'UTR, then for example if reads are 75bp, then could be the case that only 3'UTR is mapped and not coding region? I am interested in Human Sapiens. Thank you in advance.

Hi @user could you supply more detail please? A 3' UTR in viruses is a very specific region of the genome and has a precise biological meaning and purpose of study. I doubt this is your investigation, but please edit your post to be clearer.

Is it better now?Thank you!

Lets see. If you have no take up then it would be an idea to revise the question to fully explain your study. Generally it is seen as bad form to provide minimal detail, but sometimes a breif specific question connects with a brief specific answer.  This is not my area so it is difficult for me to judge.

Simple answer: Yes, that can very well happen depending on the average fragment size, the read length and the length of the UTR.


text:
Recommended coverage depth for ChIP Seq (H3K9me3 )
According to ENCODE recommendations

For broad-peak histone experiments, each replicate should have 45 million usable fragments

According to my understanding, a usable fragment means a uniquely mapped read / pair of reads.
However

H3K9me3 is an exception as it is enriched in repetitive regions of the genome. ... Tissues and primary cells should have 45 million total mapped reads per replicate.

I fail to grasp Encode's language.

Does "total mapped reads" mean uniquely mapped reads?
If so, why is 45 million "an exception" ? It is the same as the other broad peaks.

Any ideas where else can I look for an explanation?

The ENCODE Data Coordination Center, that developed and maintains the ENCODE portal, have implemented scripts validating adherence to the standards the ENCODE Consortium agreed upon when evaluating ChIP-seq analysis results. For any ChIP-seq experiment, except for H3K9me3, the adherence to the read depth standard is assessed looking on the number of "usable fragemnts" (https://www.encodeproject.org/data-standards/terms/#read-depth) (see answer by James Hawley). For H3K9me3 the read depth is assessed not by looking on usable fragments, but by looking on the total number of mapped reads (if it is a PE sequencing run, the number of mapped read pairs). It is an exception, because the number of mapped reads, especially for H3K9me3 is substantially higher than the number of usable fragments (where the duplicated and low mapping quality reads are filtered out). The decision to apply different read depth standard in case of H3K9me3 came up due to the logic explained here (https://www.encodeproject.org/chip-seq/histone/), briefly: because the H3K9me3 mark is enriched in repetitive regions of the genome, the number of usable fragments (post-filtering) is not expected to be as high as other marks (45M), despite high-quality libraries and sufficient sequencing depth. Instead of marking these experiments as failing read depth standard, the script checks the total number of mapped reads (or read pairs for PE sequencing runs) and makes sure it is above the read depth threshold (45M).

I think that "usable fragments" should also be uniquely mapped. Although the definition in Encode sates only "Fragments are considered usable if they pass the various filters in the ChIP-seq uniform processing pipelines", when I go to the filtering script at Encode's github https://github.com/ENCODE-DCC/chip-seq-pipeline/tree/master/dnanexus/filter_qc  I see "Remove unmapped reads, non-primary alignments, reads failing platform, duplicates, multi-mapping reads, and PCR duplicates."

In the initial description for "usable fragments", ENCODE has the following (emphasis mine):

Usable fragment – A fragment is defined as the sequencing output corresponding to one location in the genome. If single-ended sequencing is performed, one read is considered a fragment. If paired-ended sequencing is performed, one pair of reads is considered a fragment. Fragments are considered usable if they pass the various filters in the ChIP-seq uniform processing pipelines. Used to evaluate ChIP-seq data.

Aligning reads to repetitive regions is difficult because they don't give unique alignments.
Reads from these regions will often have poor alignment scores and may not pass alignment QC checks.
Reads may also be counted as multi-mapped, and often methods call for uniquely mapped reads, which again, may not be true for a large number of reads from H3K9me3 ChIP-seq.
So for this particular histone mark, ENCODE loosens the definition for what you should count towards a minimum sequencing depth threshold.

Why is the more difficult alignment a reason to loosen the definition? It will result in lower-quality analysis for H3K9ME3.

Btw, if you have a reference for your explanation, or can tell me how can I find one, it'd be great.

You're right, but the tradeoff seems to be between not having any usable information from H3K9me3 experiments and having a good amount of slightly lower-quality information, compared to other histone modifications.
I don't have a specific reference for that point, aside from the ENCODE documentation itself, but ENCODE has a list of all its related publications that you could search to find one.

So to be 100% sure, a "Usable fragment" is a unique alignment, right? "Corresponding to one location in the genome" means corresponding to only one location in the genome?

Yes. As per ENCODE: "Usable reads – A fragment is considered “usable” if it uniquely maps to the genome and remains after removing PCR duplicates (defined as two fragments that map to the same genomic position and have the same unique molecular identifier)."


text:
Filtering out reads from a reference (e.g. rRNA) using bowtie2
As I understand it, bowtie2 can easily be used to split reads into one of two groups:

reads for which both of a pair align well to a reference (using e.g. --al-conc-gz)
reads for which one or both of a pair do not align a reference (using e.g --un-conc-gz)

But I really want to split this second group into reads for which neither of a pair align to the reference.
bowtie2's report give a nice breakdown, but it's not clear to me how to extract specific members of these groups:
1000000 reads; of these:
  1000000 (100.00%) were paired; of these:
    958118 (95.81%) aligned concordantly 0 times
    1329 (0.13%) aligned concordantly exactly 1 time
    40553 (4.06%) aligned concordantly >1 times
    ----
    958118 pairs aligned concordantly 0 times; of these:
      179 (0.02%) aligned discordantly 1 time
    ----
    957939 pairs aligned 0 times concordantly or discordantly; of these:
      1915878 mates make up the pairs; of these:
        711518 (37.14%) aligned 0 times
        32134 (1.68%) aligned exactly 1 time
        1172226 (61.18%) aligned >1 times
64.42% overall alignment rate

In this case, I really want those that didn't align at all (i.e. the 35.58% not included in the 64.42% overall alignment rate which is the same as the 37.14% of the pairs that did not align concordantly).
Since in my case the reference is rRNA and many of those sequences are shorter than my reads, I'm not surprised that so there are so few concordant reads. I'm really interested in reads that aren't concordant and aren't discordant discordant either.
So, if I'm sticking to bowtie2, it looks like my only option is to take the resulting SAM/BAM file and filter what I want based on flags and then convert that back into FASTQ files.
Is there another option that I'm missing? It seems that the remaining output options relate to unpaired reads or SAM files.
I'm looking into other aligners now...
My goal is to go from FASTQ files to FASTQ files without extra steps, but I can go through a SAM/BAM file stage if I need to and even build a solution around that if needed.

If you want to remove rRNA, use tools that are made for this. We usually go for SortMeRNA, which can do exactly what you describe your goal is, i.e. go from FASTQ to FASTQ, splitting your reads into rRNA and non-rRNA.
In our experience, it is very accurate, if a bit slow at times.

Bam flags will indicate which reads have both pairs not mapping.

Nothing mapped.  You need to figure out why, and you can look at just about anything out figure that out.

Yes, thanks for mentioning SAM/BAM flags. I thought I mentioned them, but it looks like I only implied them. (Now updated to mention them)

I'm trying to filter out rRNA reads. So about 2/3 of the read pairs had at least one read align to my rRNA reference. However, I want the remaining 1/3 of the read pairs for which neither read aligned at all.

You are trying to filter out rRNA read by aligning them to a reference of rRNAs?  This seems ill-advised.  Filter out rRNA after aligning to your whole genome.


text:
Genome assembly of SRR12196449 with SPAdes
I am trying to assemble the run SRR12196449 with SPAdes. The description of their project is:

This project expected to standardize a method for amplification and
sequencing of FIV genome in a simple way, allowing a broader analysis
to increase the knowledge on the biology and evolution of the virus
and virus-host interaction.

This is from School of Veterinary Medicine; University of Sao Paulo.
I have tried in three ways:

Using a fasta file downloaded from NCBI (using Filtered Download method)

SPAdes-3.14.0/spades.py    --only-assembler -s sra_data.fa -o raw-fiv1

Using a fastq file also downloaded from NCBI using Filtered Download method.

SPAdes-3.14.0/spades.py   --12 sra_data.fastq.gz -o raw-fiv1-b

Downloading the original file and dumping to a fastq using sratools.

SPAdes-3.14.0/spades.py    --12 SRR12196449.fastq -o raw-fiv1-c

In the third case I get a very bad assembly, likely because it needs trimming. In the two former cases, I get an assembly with a max contig length of around 4k-5k bps. But when I use quast to evaluate the assembly I obtain a largest alignment of ~500bps. I use this as reference.
Furthermore, if make a blast search it finds that the 4k-5k contig matches other FIV sequences at ~90% . In between these other sequences, there is at least one (MF370550.1) submitted by the University of Sao Paulo (likely another run from another sample in the same project).
On the other hand, if I assemble the run using another software I obtain a contig that is an almost perfect match to my reference sequence (~99%).
I have a few questions.

Am I using SPAdes correctly? Is there a better way to use it?
Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?
What other option could explain these circumstances that I could be missing?

EDIT:
Some clarifications regarding some comments from @MaximilianPress. I can confirm that I used the same reference genome in both cases. I used that command to get the quast results:
quast-5.0.2/quast.py -r raw-fiv1/sequence.fasta SPAdes-3.14.0/raw-fiv1/contigs.fasta

The other assembler is an overlap-layout-consensus algorithm that I implemented. The methodology is similar: I use exactly the same input fasta file.
I have also tried using the SPAdes --plasmid flag and that doesn't improve much the result. Max contig length is now 3642 bps and largest alignment 737 bps, covering only 26% of the reference.
BTW, I made also tests with other virus sequences. I am just interested to know if I am using SPAdes correctly so that I can compare the performance of my algorithm Here there is a link to additional cases with other genomes. There you can find all the data and results used (raw files, assembly results and quast results)
Another clarification: I have also got perfect assemblies using SPAdes for other runs corresponding to other viruses... I am following exactly the same procedure in this case.

Can you confirm that your quast run is set up correctly? What is the other assembler tool that gives perfect results? Are there pertinent differences in methodology with the other assembler? Are there any clues in SPAdes output? Have you tried --plasmid flag? Consider reading more from the SPAdes authors about viral assembly, e.g. here: https://academic.oup.com/bioinformatics/article-abstract/36/14/4126/5837667?redirectedFrom=fulltext.

Thank you very much for the reference - that is really useful. If you could post this as an answer that would be great.

@Michael is that directed to me about posting as an answer? I gotta say, I don't feel like I answered anything!

Yes @MaximilianPress you may not feel you've answered anything ... from my point of view (which will be very similar to the OP) you definitely have! Its probably not your field (possibly)

Hey @MaximilianPress, I can confirm that I used the same reference genome in both cases. I used that command:
quast-5.0.2/quast.py -r raw-fiv1/sequence.fasta SPAdes-3.14.0/raw-fiv1/contigs.fasta
The other assembler is an overlap-layout-consensus algorithm that I implemented. The methodology is similar: I use exactly the same input fasta file. I am checking now with --plasmid flag.  Thanks for the link but it's paywalled and it seems a tutorial for identifying viruses in metagenomic samples, which is not my case... I would take a look anyway if I had not to pay.

I have tried now with --plasmid and that doesn't improve much the result. Max contig length is now 3642 bps and largest alignment  737 bps, covering only 26% of the reference.

BTW, I made also tests with other virus sequences. I am just interested to know if I am using SPAdes correctly so that I can compare the performance of my algorithm Here there is a link to more cases with other genomes link

I have also got perfect assemblies using SPAdes for other runs... I am doing the same in this case.

Hi @juanjo75es can you please update your question this this information?

Thanks @Michael, I updated it. I am not sure if that's the information you wanted to be updated...

Update 2:
I looked into this a little more, with the various data sources.
This is related in part to the answer submitted by OP juanjo75es, in addition to discussion on chat. I don't entirely understand the logic, but the general thrust seems to be that SPAdes makes weird assemblies somehow.
Some notes that I made:
REFERENCE ASSEMBLIES

FIV sequence U11820.1 was deposited in 1996, before SPAdes existed. Unclear assembly method.
FIV sequence MN630242 was deposited in 2020, but was assembled with CLC workbench.
These two sequences align together ok with minimap2. Low exact identity, it's true. This is the PAF output:

MN630242.1      8977    194     8964    +       U11820.1        9462    575     9326    1130    8778    60      tp:A:P  cm:i:116        s1:i:1125
       s2:i:55 dv:f:0.1761

These two assemblies are syntenic, and both have similar gene calls using prokka. GFFs:

MN630242.1      Prodigal:2.6    CDS     253     1605    .       +       0       ID=AMGANPBD_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00001;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     1656    4868    .       +       0       ID=AMGANPBD_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=AMGANPBD_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
MN630242.1      Prodigal:2.6    CDS     4861    5616    .       +       0       ID=AMGANPBD_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00003;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     5891    8461    .       +       0       ID=AMGANPBD_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00004;product=hypothetical protein
MN630242.1      Prodigal:2.6    CDS     8626    8790    .       +       0       ID=AMGANPBD_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=AMGANPBD_00005;product=hypothetical protein

U11820.1        Prodigal:2.6    CDS     634     1983    .       +       0       ID=EANIPDKN_00001;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00001;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     1995    5246    .       +       0       ID=EANIPDKN_00002;eC_number=3.6.1.23;Name=dut;gene=dut;inference=ab initio prediction:Prodigal:2.6,similar to AA sequence:UniProtKB:Q2YRG4;locus_tag=EANIPDKN_00002;product=Deoxyuridine 5'-triphosphate nucleotidohydrolase
U11820.1        Prodigal:2.6    CDS     5239    5994    .       +       0       ID=EANIPDKN_00003;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00003;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     6269    8830    .       +       0       ID=EANIPDKN_00004;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00004;product=hypothetical protein
U11820.1        Prodigal:2.6    CDS     8904    9152    .       +       0       ID=EANIPDKN_00005;inference=ab initio prediction:Prodigal:2.6;locus_tag=EANIPDKN_00005;product=hypothetical protein

I can share FAA files of the proteins if needed.
READ ASSEMBLY
I also assembled the indicated reads using SPAdes. For reference it is a ~9Kbp virus, but this is a 4.3Mbp assembly. There is a lot of non-virus sequence in there. The second-largest contig is a shuttle vector. the third-largest contig is cat (host). Many more are cat, so I think it's pretty cat-oriented. The original authors used CLC workbench, so I guess that just worked a lot better in this instance, even in the presence of all the contamination. Unclear why, it appears that CLC works similarly to OP's assembly tool.
I aligned these to MN630242.1 with minimap2. If I understand OP, they are unhappy about the mapping of these contigs to this genome reference. The identities are mostly high (though there is indeed some oddness with the largest contig, which only finds very low coverage- at the same time that there are overlapping contigs with very high ID?), and they cover the entire reference genome:
MN630242.1  8977    2155    3199    -   NODE_14_length_1054_cov_3786.620280 1054    2   1046    1021    1044    60  tp:A:P  cm:i:189    s1:i:1021   s2:i:378    dv:f:0.0031
MN630242.1  8977    7893    8626    -   NODE_25_length_743_cov_3814.258117  743 8   739 713 733 60  tp:A:P  cm:i:129    s1:i:713    s2:i:318    dv:f:0.0045
MN630242.1  8977    72  4904    -   NODE_1_length_4942_cov_25.814123    4942    34  4863    628 4832    41  tp:A:P  cm:i:7s1:i:628  s2:i:512    dv:f:0.1703
MN630242.1  8977    3376    3897    +   NODE_39_length_526_cov_5254.155388  526 3   524 512 521 0   tp:A:S  cm:i:8s1:i:512  dv:f:0.0023
MN630242.1  8977    6270    6673    +   NODE_82_length_409_cov_1.744681 409 3   406 357 403 0   tp:A:P  cm:i:53 s1:i:357    s2:i:351    dv:f:0.0185
MN630242.1  8977    6852    7223    -   NODE_87_length_381_cov_23.334646    381 6   376 355 371 2   tp:A:P  cm:i:5s1:i:355  s2:i:351    dv:f:0.0079
MN630242.1  8977    6187    6558    -   NODE_88_length_380_cov_2338.675889  380 7   378 351 371 0   tp:A:S  cm:i:5s1:i:351  dv:f:0.0085
MN630242.1  8977    6852    7308    +   NODE_56_length_471_cov_0.933140 471 15  471 351 456 0   tp:A:S  cm:i:43 s1:i:351    dv:f:0.0397
MN630242.1  8977    7332    7707    +   NODE_86_length_389_cov_3506.202290  389 5   380 334 375 17  tp:A:P  cm:i:5s1:i:334  s2:i:309    dv:f:0.0111
MN630242.1  8977    5181    5568    -   NODE_60_length_464_cov_1.005935 464 50  437 314 387 17  tp:A:P  cm:i:41 s1:i:314    s2:i:282    dv:f:0.0357
MN630242.1  8977    7277    7707    +   NODE_76_length_438_cov_0.919614 438 5   435 309 430 0   tp:A:S  cm:i:40 s1:i:309    dv:f:0.0437
MN630242.1  8977    6953    7285    +   NODE_94_length_345_cov_2.903670 345 8   340 301 332 0   tp:A:S  cm:i:50 s1:i:301    dv:f:0.0110
MN630242.1  8977    5744    6056    -   NODE_113_length_316_cov_2.169312    316 4   316 289 312 35  tp:A:P  cm:i:4s1:i:289  s2:i:244    dv:f:0.0096
MN630242.1  8977    5615    5927    +   NODE_108_length_322_cov_1103.400000 322 6   318 283 312 58  tp:A:P  cm:i:4s1:i:283  s2:i:209    dv:f:0.0143
MN630242.1  8977    6543    6813    -   NODE_147_length_279_cov_2592.519737 279 8   278 261 270 21  tp:A:P  cm:i:4s1:i:261  s2:i:236    dv:f:0.0032
MN630242.1  8977    8544    8872    +   NODE_100_length_333_cov_2259.189320 333 1   329 252 328 0   tp:A:P  cm:i:3s1:i:252  s2:i:245    dv:f:0.0318
MN630242.1  8977    7601    7853    -   NODE_200_length_255_cov_25.125000   255 0   252 250 252 12  tp:A:P  cm:i:3s1:i:250  s2:i:236    dv:f:0.0017
MN630242.1  8977    6648    6968    -   NODE_106_length_324_cov_2.538071    324 4   324 248 320 0   tp:A:P  cm:i:3s1:i:248  s2:i:248    dv:f:0.0303
MN630242.1  8977    5940    6193    -   NODE_162_length_266_cov_1809.287770 266 6   259 246 253 5   tp:A:P  cm:i:4s1:i:246  s2:i:240    dv:f:0.0033
MN630242.1  8977    5416    5664    -   NODE_264_length_255_cov_9.531250    255 5   253 243 248 0   tp:A:P  cm:i:4s1:i:243  s2:i:242    dv:f:0.0033
MN630242.1  8977    7221    7466    -   NODE_247_length_255_cov_12.539062   255 1   246 241 245 0   tp:A:P  cm:i:4s1:i:241  s2:i:241    dv:f:0.0032
MN630242.1  8977    4830    5071    -   NODE_1121_length_245_cov_11.711864  245 4   245 239 241 0   tp:A:P  cm:i:4s1:i:239  s2:i:235    dv:f:0.0042
MN630242.1  8977    7769    8013    +   NODE_677_length_251_cov_3.951613    251 6   249 238 244 7   tp:A:P  cm:i:4s1:i:238  s2:i:230    dv:f:0.0075
MN630242.1  8977    4901    5148    +   NODE_270_length_255_cov_8.726562    255 8   255 238 247 0   tp:A:P  cm:i:4s1:i:238  s2:i:234    dv:f:0.0073
MN630242.1  8977    8705    8948    -   NODE_382_length_254_cov_5.763780    254 7   250 234 243 15  tp:A:P  cm:i:4s1:i:234  s2:i:217    dv:f:0.0062
MN630242.1  8977    6065    6307    +   NODE_410_length_254_cov_3.488189    254 1   243 227 242 10  tp:A:P  cm:i:3s1:i:227  s2:i:210    dv:f:0.0122
MN630242.1  8977    5062    5239    -   NODE_7003_length_180_cov_9.207547   180 3   180 172 177 17  tp:A:P  cm:i:2s1:i:172  s2:i:157    dv:f:0.0066
MN630242.1  8977    11  136 +   NODE_21321_length_128_cov_1195.000000   128 2   127 125 125 3   tp:A:P  cm:i:2s1:i:125  s2:i:123    dv:f:0

Only a (relatively) few contigs align well (figure), as might be expected with heavy contamination, and of those some are overlapping:

To look into the odd behavior of the biggest contig, I ran prokka on it too and found the expected genes for its position, especially the largest gene in the virus, dut (Deoxyuridine 5'-triphosphate nucleotidohydrolase). I then took the protein sequences of the 3 assemblies and aligned them with clustalo. They're all quite similar:
CLUSTAL O(1.2.4) multiple sequence alignment

NODE1_KJHFFCBH_00001         ------MEKRPEIQIFVNGHPIKFLLDTGADITILNRKEFIIGNSIENGKQNMIGVGGGK
MN630242_AMGANPBD_00002      -------------MIFVNGYPIKFLLDTGADITILNRRDFQVKNSIENGRQNMIGVGGGK
U11820.1_EANIPDKN_00002      MGTTTTLERRLEIQIFVNGHPIKFLLDTGADITILNRKDFQIGNSIENGKQNMIGVGGGK
                                           *****:*****************::* : ******:**********

NODE1_KJHFFCBH_00001         RGTNYINVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
MN630242_AMGANPBD_00002      RGTNYTNVHLEIRDENYKTQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISDKI
U11820.1_EANIPDKN_00002      RGTNYINVHLEIRDENYRMQCIFGNVCVLEDNSLIQPLLGRDNMIKFNIRLVMAQISEKI
                             ***** ***********: **************************************:**

NODE1_KJHFFCBH_00001         PIVKVRMKDPTQGPQVKQWPLSNEKIEALTEIVERLEQEGKVKRADPNNPWNTPVFAIKK
MN630242_AMGANPBD_00002      PIVKVKMKDPNKGPQIKQWPLSNEKIEALTEIVERLEKEGKVKRADPNNPWNTPVFAIKK
U11820.1_EANIPDKN_00002      PIVKVRMRDPIQGPQVKQWPLSNEKIEALTDIVERLESEGKVKRADPNNPWNTPVFAIKK
                             *****:*:** :***:**************:******.**********************

NODE1_KJHFFCBH_00001         KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQWKKQVTVLDIGDAYFTIPLDPDYAP
MN630242_AMGANPBD_00002      KSGKWRMLIDFRELNKLTEKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
U11820.1_EANIPDKN_00002      KSGKWRMLIDFRVLNKLTDKGAEVQLGLPHPAGLQMKKQVTVLDIGDAYFTIPLDPDYAP
                             ************ *****:**************** ************************

NODE1_KJHFFCBH_00001         YTAFTLPRKNNAGPGKRYVWCSLPQGWVLSPLIYQSTLDNILQPYIKQNPELDIYQYMDD
MN630242_AMGANPBD_00002      YTAFTLPRKNNAGPGRRYIWCSLPQGWILSPLIYQSTLDNIIQPFIRQNPQLDIYQYMDD
U11820.1_EANIPDKN_00002      YTAFTLPRKNNAGPGRRYVWCSLPQGWVLSPLIYQSTLDNILQPFIKQNSELDIYQYMDD
                             ***************:**:********:*************:**:*:** :*********

NODE1_KJHFFCBH_00001         IYIGSNLSRQEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
MN630242_AMGANPBD_00002      IYIGSNLNKKEHKEKVGELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWTIQQKQL
U11820.1_EANIPDKN_00002      IYIGSNLSKKEHKQKVEELRKLLLWWGFETPEDKLQEEPPYKWMGYELHPLTWSIQQKQL
                             *******.::***:** ************************************:******

NODE1_KJHFFCBH_00001         EIPERPTLNDLQKLAGKINWASQTIPKLSIKALTHMMRGDQKLDSIREWTEEAKKEVQKA
MN630242_AMGANPBD_00002      DIPEQPTLNELQKLAGKINWASQAIPDLSIKALTNMMRGNQNLNSIREWTKEARLEVQKA
U11820.1_EANIPDKN_00002      EIPERPTLNELQKLAGKINWASQTIPDLSIKELTNMMRGDQKLDSIREWTVEAKREVQKA
                             :***:****:*************:**.**** **:****:*:*:****** **: *****

NODE1_KJHFFCBH_00001         KEAIEKQAQLYYYDPNRELYAKISLVGPHQLCYQVYHKNPEQILWYGKMNRQKKRAENTC
MN630242_AMGANPBD_00002      KKAIEGQVQLGYYDPSKELYAKLSLVGPHQISYQVYQRNPERILWYGKMSRQKKKAENTC
U11820.1_EANIPDKN_00002      KEAIEKQAQLNYYDPNRGLYAKLSLVGPHQICYQVYQKNPEHILWYGKMNRQKKKAENTC
                             *:*** *.** ****.: ****:*******:.****::***:*******.****:*****

NODE1_KJHFFCBH_00001         DIALRACYKIREESIVRIGKEPVYEIPASREAWESNIIRSPYLKASPPEVEFIHAALNIK
MN630242_AMGANPBD_00002      DIALRACYKIREESIIRIGKEPKYEIPTSREAWESNLINSPYLKAPPPEVEYIHAALNIK
U11820.1_EANIPDKN_00002      DIALRACYKIREESIIRIGKEPMYEIPASREAWESNLIRSPYLKAPPPEVEFIHAALSIK
                             ***************:****** ****:********:*.****** *****:*****.**

NODE1_KJHFFCBH_00001         RALSMVQEVPILGAETWYIDGGRRQGKAARAAYWTNTGRWQVMEIEGSNQKAEVQALLMA
MN630242_AMGANPBD_00002      RALSMIKDAPILGAETWYIDGGRKLGKAAKAAYWTDTGKWQVMELEGSNQKAEIQALLLA
U11820.1_EANIPDKN_00002      RALSMIQDAPITGAETWYIDGSRKQGKAARAAYWTDTGKWQIMEIEGSNQKAEVQALLLA
                             *****:::.** *********.*: ****:*****:**:**:**:********:****:*

NODE1_KJHFFCBH_00001         LKEGPEEMNIITDSQYILNIMNQQPDLMEGIWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
MN630242_AMGANPBD_00002      LKAGPEEMNIITDSQYVINIILQQPDMMEGIWQEVLEELEKKTAIFIDWVPGHKGIPGNE
U11820.1_EANIPDKN_00002      LKAGSEEMNIITDSQYILNIINQQPDLMEGLWQEVLEEMEKKIAIFIDWVPGHKGIPGNE
                             ** * ***********::**: ****:***:*******:*** *****************

NODE1_KJHFFCBH_00001         EVDKLCQTMMIIEGDGILNKRPEDAGYDLLAAQEIHILPGEVRIVPTRTRIMLPKGYWGL
MN630242_AMGANPBD_00002      EVDKLCQTMMIIEGDGILDKRSEDAGYDLLAAKEMHLLPGEVKVIPTGVKIMLPKGYWGL
U11820.1_EANIPDKN_00002      EVDKLCQTMMIIEGEGILDKRSEDAGYDLLAAQETHFLPGEVRIVPTKTRIMLPKGHWGL
                             **************:***:** **********:* *:*****:::** .:******:***

NODE1_KJHFFCBH_00001         IMGKSSIGNKGMDVLGGVIDEGYRGEIGVIMINLSKKSTTILEKQKVAQLIILPCKHESL
MN630242_AMGANPBD_00002      IIGKSSIGSKGLDVLGGVIDEGYRGEIGVIMINLSRKSITLLEQQKIAQLIILPCKHEVL
U11820.1_EANIPDKN_00002      IMGKSSIGSKGMDVLGGVIDEGYRGELGVIMINLTKKSITILEKQKVAQLIILPCRHESL
                             *:******.**:**************:*******::** *:**:**:********:** *

NODE1_KJHFFCBH_00001         EQGEIIMNSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNIPRIVAEEI
MN630242_AMGANPBD_00002      EQGKVIMDSERGDKGYGSTGVFSSWVDRIEEAEINHEKFHSDPQYLRTEFNLPKMVAEEI
U11820.1_EANIPDKN_00002      QQGEIQMDSERGEKGFGSTGVFSSWVDRIEEAELNHEKFHSDPQYLRTEFNLPRIVAEEI
                             :**:: *:****:**:*****************:*****************:*::*****

NODE1_KJHFFCBH_00001         KRKCPLCRIRGEQVEGKLKIGPGIWQMDCTHFNGKIIIVAIHVESGLLWAQIIPQETADC
MN630242_AMGANPBD_00002      RRKCPVCRIRGEQVGGQLKIGPGIWQMDCTHFDGKIILVAIHVESGYIWAQIISQETADC
U11820.1_EANIPDKN_00002      KRKCPLCRIRGEQVGGQLKIGPGIWQMDCTHFNGKIIIVAVHVESGFLWAQIIPQETAEC
                             :****:******** *:***************:****:**:***** :***** ****:*

NODE1_KJHFFCBH_00001         TLKAIMQLVSTHNVTEIQTDNGPNFKNQKIEGLLSYMGIKHKLGIPGNPQSQALVENANS
MN630242_AMGANPBD_00002      TVKAVLQLLSAHNVTELQTDNGPNFKNQKMEGALNYMGVKHKFGIPGNPQSQALVENVNQ
U11820.1_EANIPDKN_00002      TVKALLQLICAHNVTELQTDNGPNFKNQKMEGLLNYMGIKHKLGIPGNPQSQALVENANN
                             *:**::**:.:*****:************:** *.***:***:**************.*.

NODE1_KJHFFCBH_00001         TLKVWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYTQQESLRIQDYFSQLP
MN630242_AMGANPBD_00002      TLKAWIQKFLPETTSLENALALAVHCLNFKQRGRIGGMAPYELLAQQESLRIQEFFSKIP
U11820.1_EANIPDKN_00002      TLKAWIQKFLPETTSLDNALALALHCLNFKQRGRLGKMAPYELYIQQESLRIQDYFSQIP
                             ***.************:******:**********:* ******  ********::**::*

NODE1_KJHFFCBH_00001         QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
MN630242_AMGANPBD_00002      QKLQAQWIYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIKRVPEPCALPE
U11820.1_EANIPDKN_00002      QKLMMQWVYYKDQKDKKWKGPMRVEYWGQGSVLLKDEEKGYFLVPRRHIRRVPEPCTLPE
                             ***  **:*****************************************:******:***

NODE1_KJHFFCBH_00001         GDE
MN630242_AMGANPBD_00002      GDE
U11820.1_EANIPDKN_00002      GDE
                             ***

Overall, the level of homology is quite high at the protein level. There's maybe a little evidence that this contig's protein is closer to U11820.1 than to MN630242 (though there are also several positions where U11820.1 is the outgroup). So I'm not sure why pieces of this contig are so hard to align at the DNA level without further investigation, but honestly overall of these are looking like very similar viral sequences, as might be expected.
I also aligned the big contig to both U11820.1 and MN630242 at the DNA level in 3-way clustalo alignment it is ok. Not sure what the difference is there. There is no obvious reason to prefer one reference over the other, according to my eye. When I have clustalo output a clustering solution the big contig is essentially equidistant between the two references, so I don't buy that U11820.1 is "better" as a reference:
(
MN630242.1:0.116151
,
(
NODE1_rc:0.0504856
,
U11820.1:0.0504856
):0.0656651
)
;

As I noted in the comments, it's unfortunate that QUAST and SPAdes aren't working very well. I don't know what the deal is there. But it seems like overall the assembly is quite close to both references. I don't have any intuition for why SPAdes is doing what it's doing, or why alignment at the DNA level is acting weird (lots of artificial gaps?). But I don't think it's a pathology of SPAdes particularly, it seems like it's just something weird that happened with this dataset interacting with DNA aligners. Possibly something went strangely with read deposition?
Update:
See OP's self-answer as well (and also discussion in the comments).

Based on the (extremely informative!) follow-up updates that you have provided, I think that we can tentatively answer your questions:
The tl;dr is that you can check your reads to make sure they aren't weird and you can check your output sequence to see if it's weird.

Am I using SPAdes correctly? Is there a better way to use it?

I think that you are using SPAdes correctly, based on your experience getting ~finished viral genomes in the past with the same workflow. (THis is also my experience using SPAdes on viral genomes with Illumina PE data, is that it "just works".)
It is possible that you are using some other tools in your workflow in a suboptimal way. For example, maybe your trimming is not quite correct (it does seem based on your answer that you are already trimming). The question in my mind is then whether there is some non-intuitive thing in the reads that should still be trimmed.
I would suggest running FASTQC on your reads for this run and also your reads for the other runs that have worked, to see if something weird is going on. FASTQC directly checks for adapter sequence, weird GC composition, quality scores etc. That will give you a lot of information.

Is it possible that the sequences that blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match (not because these are correct)?

This is formally possible, of course. Our main way of knowing that you're getting the right genome is... does it look like other genomes that have been generated in the same way?
If you are really worried about this, I would suggest directly inspecting the sequence to see what the matches are. If you like, you can then blast the matches to see if they
Another approach would be to annotate the genome and make sure that it "looks like" a phage genome. E.g. it has the expected genes from FIV. Prokka is very easy to run and works ok on phage.
A slightly self-serving suggestion is to use PhageTerm to check your reads against your final genomes to ensure that they behave like phage- e.g. they have termini in the correct places, etc. (it does this with test coverage.) I am a contributor to PhageTerm, so I happen to know it a bit. Likely there are other options.

What other option could explain these circumstances that I could be missing?

As I've suggested, I think that that leaves the reads. How sure are you that they are comparable to your other runs? They look like they should be sufficient from my glance at SRA, but maybe someone bumped the MiSeq while it was running. I've given suggestions above for checking with FASTQC.

Comments are not for extended discussion; this conversation has been moved to chat.

No idea why some comments have disappeared both from here and from the chat. I find that very worrying.

@juanjo75es I believe that Devon Ryan used his StackExchange superpowers to move our discussion in comments to a chat. I think that we can still see all the comments by following the link they posted.

Regardless of superpowers its a great discussion and really interesting

After many considerations, I am going to accept the response from Maximilian Press. I see now that some viruses have high variability (HIV even 50% of the sequence). Therefore MN630242.1. and U11820.1 are apparently two strains. There are things I still don't understand but these are beyond the initial goal of my question. Concretely:

Why SPAdes returns one strain and rnaSPAdes the other one.
Why one strain matches 99,9% with MN630242.1 (and at least another assembly) while the other one has the (now) expected variability of this kind of virus.

I also want to point that apparently, Quast is not effective for viruses with such high variation.
I'm not going to delete this answer given that it responds to part of my question that is not responded in Maximilian's answer.
Therefore, directly answering my questions:

Am I using SPAdes correctly? Is there a better way to use it?
I am using SPAdes mostly in the correct way, similarly to how an average user would. But rnaSPAdes seems to be more appropriate for RNA viruses and it indeed works far better in this case.

Is it possible that the sequences that Blast finds matching with the SPAdes contigs are indeed sequences that someone got using SPAdes and that's why these match?
That could be possible in some cases but that doesn't mean these sequences are wrong as I initially considered. For whatever reason SPAdes seems to miss one strain (if that's what really is happening)

What other option could explain these circumstances that I could be missing?
See Maximilian Press answer.

Uhm, it would have been helpful having known from the beginning that I could add images...

I want also to note that I don't agree with the argument about contamination in @MaximilianPress answer and that there is still an issue with all that which might be altering some studies: they could be obtaining the sequence of a random quasispecies of the virus. Seems a complex issue far from my actual interests

I think that this is the more effective answer than mine from my perspective, as it highlights the usefulness of rnaSPAdes in this context (retroviruses). My answer might be considered more useful from the general perspective of evaluating whether assemblies are plausible (arising from discussion with @juanjo75es), but that is less relevant to the practical matter of obtaining a better FIV assembly.


text:
Too slow issue of BioMart
I am using library(biomaRt) in R to annotate SNP rsid. It takes too much time for only 1000 SNPs....... I would like to know how I can run this code faster.
 coords<-data
 coords$chr_name<-coords$chromosome
 coords$start<-coords$position
 coords$end<-coords$position
 coords<-coords%>%dplyr::select(chr_name,start,end)
 mart=useMart(biomart="ENSEMBL_MART_SNP", host="grch37.ensembl.org",
 path="/biomart/martservice", dataset="hsapiens_snp")
 results=data.frame(matrix(ncol=4,nrow=0))
 for (i in seq (1,nrow(coords)))
{
ens=getBM(attributes = c('refsnp_id','allele','chrom_start','chrom_strand'),
filters = c('chr_name','start','end'), values = as.list(coords[i,]), mart = mart)
results=rbind(results,ens)}

I am using the following code, but still not working......
coords$newposition<-paste0(coords$chromosome,":",coords$position)
coords$chromosomal_region<-paste0(coords$newposition,":",coords$position)
coords<-coords$chromosomal_region
res<-getBM(attributes = c('refsnp_id','allele','chrom_start','chrom_strand'),
   filters = 'chromosomal_region',
   values = coords,
   mart = mart)

I am also facing the same issue. I have 5600 SNPs and even after vectorization, I am facing the same timeout problem even for 500 SNPs.

Don't run getBM 1000 times.  It can handle whole vectors, and get you all your results at once.

Sorry, could give me some example codes?

vectorising the function is a good idea and would speed up the process. It would be nice if @swbarnes2 gave a less terse answer with e.g some example code on how to achieve it.


text:
I want to display color key on my heatmap
library(gplots)

data=structure(c(-0.1152, 0.1548, 0.0037, 0, 0, 0, 0, 0, 0, 0.1369,
0.1914, 0, 0.0499, 0, 0, 0, 0, 0.1259, 0, 0.1854, -0.4412, 0,
-0.4425, 0, 0, -0.5141, 0, -0.3463, 0.1529, -0.1176, 0.4153,
-0.4412, 0.2169, 0.1894, 0.3288, 0.5549, -0.2549, 0.1915, -0.0759,
0.1255, 0.2269, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.4412, -0.1923,
0.1502, 0, 0, 0, -0.1978, 0.0828, 0, 0.2214, 0.192, 0, -0.2651,
0.192, -0.1412, -0.287, 0, -0.2299, 0, 0.1095, 0, 0.0166, 0,
0.1617, -0.16, 0, 0, 0.1255, 0, 0.2298, 0.153, 0.1627, 0, 0.293,
0.2476, -0.1565, 0.0899, 0, 0.2217, 0.2274, 0.1764, 0.1538, 0,
0.1747, 0.1857, -0.2309, 0.0827, -0.0609, 0.0754, -0.2653, 0.1461,
0.1611, -0.1201, 0, 0, 0.114, 0.237, 0.1738, 0, 0.2269, 0, 0.0068,
0, 0.1288, -0.1727, -0.1187, -0.1502, -0.1502, 0, -0.4104, -0.1305,
-0.1126, -0.2642, -0.1264, -0.2403, 0.1522, 0.0211, 0.1175, -0.0549,
-0.2571, -0.1305, -0.2143, -0.2153, -0.0566, -0.2068, 0.0406,
-0.087, 0.079, -0.0238), .Dim = c(10L, 14L), .Dimnames = list(
    c("AI", "AUDC", "AUDT", "AUT", "BD", "CB", "CD", "DPW", "MDD",
    "SZ"), c("BI", "BB", "CD", "CGC", "CTF", "HDT", "MTS", "NOI",
    "OFI", "SRV", "SAF", "FFI", "SUM", "WIN")))

data <- read.csv("C:/Users/Lenovo/Desktop/Correlation matrix VITD _Disorder_R/Z1.txt", header=TRUE, row.names=1,sep="")
data = as.matrix(data[![enter image description here][1]][1])
heatmap.2(data,
          trace ="none",
          density = "none",
          col = redblue(100),
         key=TRUE,
         keysize=1,
          cexRow =1,
          cexCol =1,
          margins = c(5,5),
          labCol = colnames(data),
          scale="row",
          hclust=function(x) hclust(x,method = "complete"))

  [1]: https://i.stack.imgur.com/jlwf4.png

I ran your code and it already displays a colour key - what's the issue?]

Ah @user438383 your edit makes the data look alot better. I suspect that will now run fine. You could just post your edit as an answer. In fact I would  recommend it because otherwise a resolved question will be scored as an unresolved, unanswered question.

I'm not sure exactly what was going on with your data, but I read it in as a raw data.frame using as.matrix(read.table(., row.names=1, header=T)). The issue was that you were reading it in without a header, so the col.names were being read in as part of the matrix. But apart from that, the rest of your command was fine:
library(gplots)

data = structure(c(-0.1152, 0.1548, 0.0037, 0, 0, 0, 0, 0, 0, 0.1369,
0.1914, 0, 0.0499, 0, 0, 0, 0, 0.1259, 0, 0.1854, -0.4412, 0,
-0.4425, 0, 0, -0.5141, 0, -0.3463, 0.1529, -0.1176, 0.4153,
-0.4412, 0.2169, 0.1894, 0.3288, 0.5549, -0.2549, 0.1915, -0.0759,
0.1255, 0.2269, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -0.4412, -0.1923,
0.1502, 0, 0, 0, -0.1978, 0.0828, 0, 0.2214, 0.192, 0, -0.2651,
0.192, -0.1412, -0.287, 0, -0.2299, 0, 0.1095, 0, 0.0166, 0,
0.1617, -0.16, 0, 0, 0.1255, 0, 0.2298, 0.153, 0.1627, 0, 0.293,
0.2476, -0.1565, 0.0899, 0, 0.2217, 0.2274, 0.1764, 0.1538, 0,
0.1747, 0.1857, -0.2309, 0.0827, -0.0609, 0.0754, -0.2653, 0.1461,
0.1611, -0.1201, 0, 0, 0.114, 0.237, 0.1738, 0, 0.2269, 0, 0.0068,
0, 0.1288, -0.1727, -0.1187, -0.1502, -0.1502, 0, -0.4104, -0.1305,
-0.1126, -0.2642, -0.1264, -0.2403, 0.1522, 0.0211, 0.1175, -0.0549,
-0.2571, -0.1305, -0.2143, -0.2153, -0.0566, -0.2068, 0.0406,
-0.087, 0.079, -0.0238), .Dim = c(10L, 14L), .Dimnames = list(
    c("AI", "AUDC", "AUDT", "AUT", "BD", "CB", "CD", "DPW", "MDD",
    "SZ"), c("BI", "BB", "CD", "CGC", "CTF", "HDT", "MTS", "NOI",
    "OFI", "SRV", "SAF", "FFI", "SUM", "WIN")))

heatmap.2(data,
          trace ="none",
          density = "none",
          col = redblue(100),
         key=TRUE,
         keysize=1,
          cexRow =1,
          cexCol =1,
          margins = c(5,5),
          labCol = colnames(data),
          scale="row",
          hclust=function(x) hclust(x,method = "complete"))

Super answer, thank you and thanks for editing and clarifying the original question.

Thanks for the help. I have run the code exactly as you have sent it, but still no key.It is really annoying.

@AranyakGoswami Could you edit your main post to show the exact command you have used and post a picture of the figure you created (without the Key), as well as the results of sessionInfo()? Thanks.

Hi @AranyakGoswami this is a really good answer. Could please mark this resolved and start a new question? I wasn't familiar with the gplots package. Sectioning a project into a series of questions is the best way to resolve it and fairest way for the respondee

the issue is resolved.

Hi @AranyakGoswami Great! Please could you sign the answer as 'resolved' and a green tick will appear besides the answer? This avoids the question (which has been well answered) being recycled at a later date and gives credit accordingly.

resolved. However I am still unable to have the color key on my heatmap.


text:
error in random forest analysis
I am now struggling to do random forest analysis, I will be thankful if you could help with code for random forest analysis.
I got samples from the root, soil, and leaf from two regions (bau & mau) and these samples belong to two seasons (Wet and Dry).
Now I am interested to do random forest analysis at genera or family level to identify the taxa which contribute the differences like in root samples based on region as well as season.
Here is my code, but I am getting the error.
    library(randomForest)
    library(knitr)

    #### RANDOM FOREST ANALYSIS #####
     #### Prepare data ####
     #Load OTU table
     OTU_table=t(read.table("asv.table.txt", row.names=1,sep="\t", header=T, blank.lines.skip=F, check.names=F))
     table(apply(OTU_table,1,sum)) #verify rarefaction

     #Load metadata
     Meta=read.table("metadata.txt", header=T, row.names=1, stringsAsFactors=F, na.strings="NA",check.names=FALSE)
     Meta$sampleid=rownames(Meta)

     #Load taxonomy
     taxo=read.table("taxonomy.txt", row.names=1, sep="\t", header=F ,stringsAsFactors=F,quote="")
     rownames(taxo)=paste("a.",row.names(taxo),sep="")

     #### Run models ####
     #1. Root only
     # 1.1. both region
     # 1.2. bau only
     # 1.3. mau only
     #2. Soil only
    # 2.1. both region
     # 2.2. bau only
     # 2.3. mau only
     #Params RF
     NTREE=1000 # Number of Trees
     NbVar=1000 # Number of variables tested at each split

     #### Root ONLY 1-3 ####

 # 1. BOTH Region

 #Subset of data
 RootSamples=as.character(Meta[Meta$Compartment=="Root","sampleid"])
 Root_OTU_table=OTU_table[RootSamples,]

 #Model with microbiome based on Season, region

whole_root_pred=data.frame(Season=Meta[RootSamples,"Season"],Region=Meta[RootSamples,"Region"],a=Root_OTU_table)

head(whole_root_pred)
     Season Region a.d2ec9f3b77975c0f457e4b7413b217ff
     a.3147790f0d5a78316fb9dd64f53b9473 a.97aecc1f35cc1f50db507ad71dd22367
     a.bfad6370d28182cc6304844e9bec7fb6 a.5fa2a987221a1d9ca416148570c18086

    **RF_model_Root_all=randomForest(y=?,sampsize=c(143,143),strata=?,x=whole_Root_pred,importance = T,proximity = T,ntree =
    NTREE,mtry = NbVar)**

print(RF_model_Root_all)
     #plot summary using the 5% most important OTUs ERROR ON LAST LINE
     imp=data.frame(importance(RF_model_Root_all))
     imp$genus=as.character(taxo[rownames(imp),"Genus"])
     Best=imp[imp$MeanDecreaseAccuracy>quantile(x = imp$MeanDecreaseAccuracy,.95),]
     bymedian <- with(Best, reorder(genus, -MeanDecreaseAccuracy, median))

     pdf(width = 20,height = 10,file=paste(pathforplots,"Variable_Importance_Root_BothRegion_raref.pdf",sep=""))
     par(mar=c(15,5,1,1))
     boxplot(Best$MeanDecreaseAccuracy ~ bymedian, data = Best,
     xlab = "", ylab = "Variable Importance",
     main = paste("Root in Both Countries; Error Rate=",round(RF_model_Feces_all$err.rate[NTREE,"OOB"],3),sep=""), varwidth = TRUE,
     col = "lightgray",las=2)

     dev.off()

Many thanks

What was the outcome for the question on alpha diversity,
https://bioinformatics.stackexchange.com/questions/14501/alpha-diversity-wilcox-test ?

This question will be a little hard to answer without more information.
For example, we will need to see your dataset (whole_root_pred), to decide why Stunting_Root is NULL.

You might need to initialize Stunting_Root as a variable. It is currently not clear if it is e.g. a column of your dataframe, or just uninitialized. Uninitialized variables are NULL, which would explain your problem. randomForest might not know to look for strata inside your dataframe, for example. Is it in your dataframe?
Also, I might be missing something, but why are you passing ? as a response? I'm not an expert but I believe that is an illegal character in R (I'm pretty sure?).

Thank you so much for helping.  I have updated the code. Actually after subset root data, I am not sure what to put in RF modal command in Y=? , as well as in strata?. Till whole_root_pred command, I am not getting any error but not sure after that what to do?

@bioinfonext I would recommend looking at the randomForest docs (?randomForest) to decide what to put in these arguments. y is supposed to be the response variable, so either a vector of categories (if task is classification) or a vector of numerics (if task is regression). Again, it would be very helpful to see a sample of the data, both the response variable and the predictor variables, and how they are represented in the code you provide.

Thanks how can I show sample data! Is it head or dput command?

head() should be fine.

I have added head(whole_root_pred) output.


text:
project ideas for a final year CS undergrad/CS masters student
I am enrolled in an online undergrad CS program in Canada and also Georgia Tech OMSCS where I finished computational photography and AI for Robotics.  I also have some deep learning programming knowledge.
For my capstone project in the undergrad, I was building a nutrition app but I lost interest after finding so many other apps and repos with similar ideas.  So I though since I'm taking an intro bioinformatics course this final term before I graduate, I could do a bioinformatics programming project.
Does anyone have any good ideas?  I'm free labor if I can get approval from my course supervisor for your proposal.  Time length would be approximately 10 weeks.
My languages in order of familiarity: Python, C++, C#, Java, C, Javascript, HTML, CSS.

I’m voting to close this question because it is off-topic and unrelated to bioinformatics. In general you should talk to experienced supervisors to guide you for project ideas that can be finished it time while neither over- nor underdemanding you.


text:
panic: division of zero by zero or infinity by infinity error when running GODON positive selection analysis
I am using godon to find genes under positive selection with the branch-site test. However, after running for a little while I get the error "panic: division of zero by zero or infinity by infinity".
Below is the full output with the error report. I can supply the phylogenetic tree and fasta file I am using if that would be helpful. Thank you.

iter=33 lnL=-73104.072
iter=66 lnL=-72956.855
iter=99 lnL=-72912.956
iter=131 lnL=-72902.900
iter=162 lnL=-72900.147
iter=194 lnL=-72898.562
iter=226 lnL=-72897.970
iter=258 lnL=-72897.651
iter=289 lnL=-72897.435
iter=321 lnL=-72897.284
iter=353 lnL=-72897.199
iter=385 lnL=-72897.153
iter=416 lnL=-72897.108
iter=448 lnL=-72897.076
iter=17 lnL=-72890.551
panic: division of zero by zero or infinity by infinity
goroutine 1 [running]:
math/big.(*Float).Quo(0xc4273a8f40, 0xc4273a8f68, 0xc4273a8f18, 0xc42e88c118)
/usr/local/go/src/math/big/float.go:1639 +0x15c
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).BEBPosterior(0xc422164120, 0xc4259bb900, 0x91, 0x91)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:416 +0x6e9
bitbucket.org/Davydov/godon/cmodel.(*BranchSite).Final(0xc422164120, 0x101)
/go/src/bitbucket.org/Davydov/godon/cmodel/bs.go:455 +0x112
main.computeFinal(0x1f263a0, 0xc422164120, 0xc4201522a0, 0x101, 0x6a80, 0x0)
/go/src/bitbucket.org/Davydov/godon/godon/optimize.go:150 +0x93
main.performSingleTest(0xc420162320, 0xc428266000, 0x9399, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, 0x0, ...)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:335 +0x1b52
main.hypTest(0x1ee03a7, 0x1ee03a7, 0x4, 0xc4201aa501, 0x1, 0x1)
/go/src/bitbucket.org/Davydov/godon/godon/hyptest.go:88 +0x71a
main.main()
/go/src/bitbucket.org/Davydov/godon/godon/godon.go:316 +0xefa

Hi, sorry about delay. Apparently I didn't subscribe properly to the tag.  Anyhow this looks like a bug to me, would be great if you could open an issue on bitbucket and attach an example file together with options you used to run software. https://bitbucket.org/Davydov/godon/issues

Ah wait you mentioned a branch sites test, right that makes sense. You've got two or more identical sequences in the analysis. Ensure there are no identical sequences in your alignments and ... PRESTO! It will work :-)
Below are my original thoughts before I re-read your question. I've kept the notes because it helps explain the theory a bit. There is no need to supply the tree, it should be resolved.

The tree likelihood looks to be decreasing nicely, so it can see both the tree and the codons in the alignment.
At a guess the error relates to two or more indentical nucleotide sequences in your alignment. Thus zero non-synonymous sites (amino acid change) and zero synonymous (silent) sites. This error would only occur if you were calculating selection for each 'node' in the tree (internal branch).  Alternatively, possibly if the number of synonymous sites was zero, regardless of the non-synonymous mutations this error could appear, but this doesn't seem to fit the error output.
I don't know "Godon", but recommend PAML for this analysis. Admittedly it is complicated to deploy, but it is a very robust algorithm and hugely respected.