cannin/biostars_qa huggingface dataset:
Question: Hello All,

Most of the identified bacterial genes are hypothetical. I want to build a database of any and all bacterial genes so I can do an RNA-seq simulation of a microbiome.

However, I have yet to find a pre-build database like this. Refseq database does include several RNA.fna files, but the sequences are a small subset of what is available since it is curated.

Can anyone advise on the best way to make an query to a repository like GenBank similar to "all prokaryote, nucleotide, protein encoding" sequences.

Any solution is good, even an un-elegant one, such as downloading and parsing through all the database files myself. I could really use guidance on how to approach the problem.

Answer: I'm not sure of your rationale for creating your own RNA-seq simulation when there are numerous mock and real data metatranscriptomes for many environments, but especially for human microbiome (you didn't clarify -- the expressed reads would be different for differing microbiomes).

If you're looking for a master list of all bacterial coding regions, I would just pull all the bacterial genomes from NCBI using their FTP server and select all the protein encoding genes from the genomes using the GFF files. This shouldn't take long and probably downloading all the data would be the most time consuming part. Selectively downloading just the coding regions would be the best use of your time -- not sure if NCBI has a list of files. Searching FTP is troublesome.
Question: Hi everyone,

I am using the [EstimateLibraryComplexity][1] utility from Picard Tools to calculate the complexity of my paired-end RNAseq libraries.

[1]: http://picard.sourceforge.net/command-line-overview.shtml#EstimateLibraryComplexity

This is my command line:

java -jar /picard/EstimateLibraryComplexity INPUT=sample.bam OUTPUT=sample_libcomp.txt VERBOSITY=ERROR VALIDATION_STRINGENCY=SILENT

This generates a `sample_libcomp.txt` file. This is the truncated output:

```
## HISTOGRAM java.lang.Integer
duplication_group_count P01311

1 23739815
2 3633946
3 870509
4 426481
5 202751
6 171461
7 93221
8 83632
9 58171
10 50066
11 34938
12 36788
13 24277
14 24100
15 19388
16 18345
17 13640
18 14480
...
456 1
457 1
458 1
459 1
460 2
464 3
468 1
470 2
471 2
473 1
477 2
480 1
484 1
488 1
```

Can anyone explain to me what these values mean? I couldn't find an explanation of the output anywhere. I plan to plot these values as a density histogram (maybe convert the values to log2). So I really need to understand what these values are in order to interpret the histogram that I will create later.

Thanks!

Answer: The first column is the number of duplicates. The second column is the number of reads having the corresponding number of duplicates.

So in your output, there are 426,481 sequences which have exactly four duplicates.
Question: I'm trying to run remote blast, but I keep getting this error:

**-bash-4.2$ blastp -remote -db swissprot -query test.txt -out outfile.txt -outfmt 7

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (310.5) [blast4] Service not found

Error: (315.2) CConn_Streambuf::CConn_Streambuf(): NULL connector: Unknown

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (310.5) [blast4] Service not found

Error: (315.2) CConn_Streambuf::CConn_Streambuf(): NULL connector: Unknown

Error: (802.5) CObjectOStream: error at byte 14: stack is empty: write fault

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (310.5) [blast4] Service not found

Error: (315.2) CConn_Streambuf::CConn_Streambuf(): NULL connector: Unknown

Error: (802.5) CObjectOStream: error at byte 14: stack is empty: write fault

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (303.7) [URL_Connect] Failed to connect to www.ncbi.nlm.nih.gov:443: Not supported

Error: (310.5) [blast4] Service not found

Error: (315.2) CConn_Streambuf::CConn_Streambuf(): NULL connector: Unknown

Error: (802.5) CObjectOStream: error at byte 0: stack is empty: cannot close output stream

Error: (CIOException::eFlush) COStreamBuffer::Flush: failed

Error: (CSerialException::eIoError) byte 0: cannot close output stream

Error: (802.4) Cannot close serializing output stream(CSerialException::eIoError) byte 0: cannot close output stream

Error: (802.5) CObjectOStream: error at byte 14: stack is empty: write fault

Error: NCBI C++ Exception:

"/opt/insy/env.el7/src/ncbi-blast-2.2.27+-src/c++/src/util/strbuffer.cpp", line 794: Error: ncbi::COStreamBuffer::FlushBuffer() - write fault

"/opt/insy/env.el7/src/ncbi-blast-2.2.27+-src/c++/src/serial/objostr.cpp", line 571: Error: ncbi::CObjectOStream::Write() - Frame type= eFrameNamed, Object type= Blast4-request*

I am pretty sure that my command is correct and I have a working internet connection. Any ideas?
Thanks!

Answer: Have you downloaded the latest blast+ package (v.2.5.0)? If not that should be the first thing to try.

NCBI has moved to `https-only` connections, which are only supported in the latest versions of their software (e.g. blast, sratoolkit etc). You appear to be using blast+ v.2.2.27.
Question: Hi,

Section "2.10 Record inputs" in the UserGuide is a bit unclear to me. "dependent_parameters" is used as an id and as a name. Which is the actual keyword and which is an arbitrary name that I can change ?

The example shows how to make one set of dependent parameters. How can I make more than one set ? E.g.: itemA pairs with itemB and itemD pairs with itemE.

thanks
Wolfgang

Answer: Hello Wolfgang,

Neither usage of `dependent_parameters` is a keyword, in fact there is an error in the v1.0 spec related to custom records and the `name` field: https://github.com/common-workflow-language/common-workflow-language/issues/394

You can have as many dependent parameters as you want, just follow the pattern: custom record type with two or more fields.

Here's a real example that may help:

External definition of custom record type for re-use:

type: record
name: model
fields:
- name: main
type: File
- name: prob_forward
type: File
label: Forward probability model
[…]
(From https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/f8bb388/tools/FragGeneScan-model.yaml)

Import of the custom record type in the tool definition

requirements:
SchemaDefRequirement:
types:
- $import: FragGeneScan-model.yaml
[…]
(From https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/f8bb388/tools/FragGeneScan1_20.cwl#L21)

Use of the custom record type in that tool's inputs:

inputs:
[…]
model: FragGeneScan-model.yaml#model
(From https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/f8bb388/tools/FragGeneScan1_20.cwl#L51)

And likewise, using that tool in a workflow. First import the same custom type:

requirements:
[…]
- class: SchemaDefRequirement
types:
- $import: ../tools/FragGeneScan-model.yaml
# other imports for other tool's custom types go here
(from https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/f8bb388/workflows/emg-pipeline-v3.cwl#L9)

As needed reference it -- for example in the workflow level inputs:

inputs:
[…]
fraggenescan_model: ../tools/FragGeneScan-model.yaml#model
(From https://github.com/ProteinsWebTeam/ebi-metagenomics-cwl/blob/f8bb388/workflows/emg-pipeline-v3.cwl#L18)

Thank you for your question, this part of the user guide is confusing and needs rewording. There is an open issue to remind the team to clarify and demonstrate this better: https://github.com/common-workflow-language/common-workflow-language/issues/438
Question: Hi :

I have a file of Data includes 20 proteins, each 10 proteins belongs to different class. My goal is to extract the similarities among 10 sequences and between those two classes in Protein level and RNA level. So, I reversed them to RNA code and I saved them in different file. I want to analysis both files of protein and RNA to find some similarity regions that may RNA share it in each class, and repeat the same thing in Protein. Firstly, I use local multiple sequence alignment of 10 RNA sequences by MUSCLE , and I used Jalview program for this purpose. Jalviews shows me some coloured area that have the same percentage of identity.

My question is:

I would like to represent these colours by numbers, but Jalview didn't give me any score or numbers to find the percentage of identity.. !! How I can analysis my data to extract the similarities among sequences of both class and both level of gene (RNA and Protein) ? And how I can represent them by quantitative measurement (e.g. similarity score, level of variation, percentage of identity) or you can suggest me another type of measurement ?

Thanks ..

Answer: Consurf may do what you want:

http://consurf.tau.ac.il/overview.html

d. Calculating the amino acid conservation scores
The conservation score at a site corresponds to the site's evolutionary rate. The rate of evolution is not constant among amino (nucleic) acid sites: some positions evolve slowly and are commonly referred to as "conserved", while others evolve rapidly and are referred to as "variable". The rate variations correspond to different levels of purifying selection acting on these sites. The purifying selection can be the result of geometrical constraints on the folding of the protein into its 3D structure, constraints at amino acid sites involved in enzymatic activity or in ligand binding or, alternatively, at amino acid sites that take part in protein-protein interactions.

In ConSurf, the rate of evolution at each site is calculated using either the empirical Bayesian [11] or the Maximum Likelihood [12] paradigm. In both of these methods, the stochastic process underlying the sequence evolution and the phylogenetic tree are explicitly taken into account. The Bayesian method was shown to significantly improve the accuracy of conservation scores estimations over the Maximum Likelihood method, in particular when a small number of sequences are used for the calculations [11] . An additional advantage of the Bayesian method is that a confidence interval is assigned to each of the inferred evolutionary conservation score.
Question: I am looking for a tool/script/pony to correct the REF column in a vcf file whenever that nucleotide doesn't match the reference genome, as supplied in fasta format.

It sounds like a common task but I could not find something. I did find `bcftools +fixref` but that only works for SNPs. My vcf files are from structural variants.

Cheers,
W

Answer: Some code like this can help, but be aware that changing the reference column affects other fields in your VCF, such as the GT value.

#!/usr/bin/perl

use strict;
use warnings;

$ARGV[2] or die "use fixRefVCF.pl \n";

my $genome = shift @ARGV;
my $orig_vcf = shift @ARGV;
my $new_vcf = shift @ARGV;

my %seq = ();
my ($chr, $pos, $ref, $obs);

warn "loading genome sequence from $genome\n";
open (my $fh, "<", $genome) or die "cannot read $genome\n";
while (<$fh>) {
chomp;
if (/>(\w+)/) {
$chr = $1;
} else {
$seq{$chr} .= $_;
}
}
close $fh;

open (my $vh, "<", $orig_vcf) or die "cannot read $orig_vcf\n";
open (my $oh, ">", $new_vcf) or die "cannot write $new_vcf\n";
while (<$vh>) {
if (/^#/) {
print $oh $_;
} else {
chomp;
my @var = split (/\t/, $_);
$chr = $var[0];
$pos = $var[1]; # VCF is 1-based, Perl string is 0-based
$ref = $var[3];
$obs = substr($seq{$chr}, $pos - 1, length $ref);
if ($ref ne $obs) {
warn "editing $chr:$pos $ref->$obs\n";
$var[3] = $obs;
}
print $oh join "\t", @var;
print $oh "\n";
}
}
close $vh;
close $oh;
Question: I have a file like this:

"

DGM97JN1_135:2:1101:1283:2110 16 chr13.fa ...

DGM97JN1_135:2:1101:1434:2186 16 chr08.fa ...

DGM97JN1_135:2:1101:1385:2244 0 chr16.fa ...

DGM97JN1_135:2:1101:1663:2038 0 chr13.fa ...

...........

"

and I would like to use a fast and easy way to delete the ".fa" after chromosome names in each row of the third column. How I do it now is

cat <input file> | tr -d '.fa'

But it does not works well, since it will also delete all "." in the file for some reasons, why this would happen and what is the right way to code it?

Besides since it is a large file, I wonder if there is a way I can narrow down the searching and substituting on only third column and therefore accelerate the process?

Answer: you need sed, not tr

try:

sed 's/\.fa//' <input file>
Question: How to choose between signed vs unsigned network ?

sft=pickSoftThreshold(datExpr,dataIsExpr = TRUE,powerVector = powers,corFnc = cor,corOptions = list(use = 'p'),networkType = "signed")

adj= adjacency(datExpr,type = "unsigned", power = softPower)

TOM=TOMsimilarityFromExpr(datExpr,networkType = "signed", TOMType = "unsigned", power = softPower)

These are the three lines where i have to defined signed or unsigned type i read through Peter Langfelder blog but still im not much clear about

Any suggestion would be highly appreciated

Answer: I think the "sign" represents the sign of weight on the edges. If you care about the sign, for example sign represents positive or negative regulation between two nodes, then you need to choose "signed". On the other hand, if the weight just represent the strength of relatedness between two nodes, you might don't care positive or negative, so you need to choose "unsigned". That really depends on your needs.
Question: Hello all,

I am looking at some mass-spec data. I found several fragments mapping to **Ig heavy chain V-II region WAH protein** and want to find corresponding gene.

Example http://www.uniprot.org/uniprot/P01770

Uniprot says the gene name as "NULL". Is this an annotation error or any special aspect of Ig regions am missing? I want to map several proteins with these type of names to genes.

- Cluster of Ig heavy chain V-I region HG3
- Cluster of Ig heavy chain V-II region SESS
- Cluster of Ig heavy chain V-III region BRO
- Cluster of Ig lambda chain V-I region NEW
- Cluster of Ig lambda chain V-II region BUR
- Ig heavy chain V-II region WAH
- Ig heavy chain V-III region BUT
- Ig heavy chain V-III region GAL
- Ig heavy chain V-III region NIE
- Ig heavy chain V-III region WEA
- Ig kappa chain V-I region Kue
- Ig kappa chain V-I region Wes
- Ig kappa chain V-III region VG (Fragment)
- Ig lambda chain V-III region LOI
- Ig lambda chain V-III region SH
- Ig lambda chain V-V region DEL

How can I map these to corresponding gene names?

Thoughts?

Answer: From the UniProt helpdesk (please contact us directly with future questions):

There seem to be 281 human UniProtKB/Swiss-Prot entries (http://www.uniprot.org/uniprot/?query=+NOT+gene%3A*+AND+organism%3A%22Homo+sapiens+%28Human%29+[9606]%22&sort=score) that do not have any gene name. This figure represents less than 1.4% of the total, which is not too bad. Sequences that have been exclusively produced by large scale sequencing programs, such as NEDO, for instance, usually do not have any satisfactory name, but rather a clone name, such as FLJxxx in the case of the NEDO project. We do not consider these names as eligible for gene names and there are no publications that could give us a hint to name these sequences properly.

There are about 65 entries of the FLJxxx type in UniProtKB/Swiss-Prot.

Although I didn't look at all entries that do not have any gene name, I expect that most follow the same rationale. At some point, these entries will be either "upgraded" if some additional evidences are produced by the scientific community, or deleted from the database. However, UniProtKB/Swiss-Prot is quite conservative and we prefer to keep some dubious sequences. This provides MS users with the largest possible set of peptides for their identifications.

There are also human entries that do not have any official gene name. Currently 356 human entries do not have any link to HGNC. This represents less than 2% of the total number of human entries. Some may belong to the group of orphan sequences that do not seem to interest anyone for the time being (see above), or they may be missing from HGNC (see for instance Ovostatin homologs, Q6IE37 and Q6IE36).
Question: I have a question concerning the extraction of sequences from a fasta file (>7000 sequences) using a reference .txt file with sequence headers. I have been playing around and been looking all over the internet to find a solution for this problem, but surprisingly, nothing really matches what I want to do. So, I have two files:

**1) a fasta file which looks like this:**

>Zotu1

ACTGACAAAGCA

TGCACGTCATTTT

>Zotu2

ATGCATCAGCATA

TGACCCCCGTTTA

>Zotu10

CGTCGAAAAATTT

CGATACACCCTAT

>Zotu22

CGTACGTCCCCTT

CGATATAATATATA

2) a .txt file with a list of sequence names:

Zotu1

Zotu2

Now, I want to use the .txt file to select sequences from the .fasta file. I have two semi-solutions that do part of the job.

**OPTION 1:**

cat list.txt | awk '{gsub("_","\\_",$0);$0="(?s)^>"$0".*?(?=\\n(\\z|>))"}1' | pcregrep -oM -f - sequences.fas > newfile.fas

*Problem:* this function gives me the full sequences, but extracts too many sequences since everything that partially matches the strings in the .txt file will be selected. In this case, it means that also Zotu10 and Zotu22 are selected.

**OPTION 2:**

grep -A 1 -wFf list.txt sequences.fas > newfile2.fas

*Problem:* this function correctly selects only the sequences that completely match the strings in the .txt file, but does not return the full fasta sequences, but only the part of the sequence on the first line. An output thus looks like this:

>Zotu1

ACTGACAAAGCA

>Zotu2

ATGCATCAGCATA

I tried combining both solutions but that somehow did not end well. I would be much helped by an elegant solution for this problem, preferably using the codes I already obtained.

Many thanks!

Answer: with seqtk:

$ seqtk subseq test.fa test.txt
>Zotu1
ACTGACAAAGCATGCACGTCATTTT
>Zotu2
ATGCATCAGCATATGACCCCCGTTTA

with grep:

$ grep -w -A 2 -f test.txt test.fa --no-group-separator
>Zotu1
ACTGACAAAGCA
TGCACGTCATTTT
>Zotu2
ATGCATCAGCATA
TGACCCCCGTTTA

inputs:

$ cat test.fa
>Zotu1
ACTGACAAAGCA
TGCACGTCATTTT
>Zotu2
ATGCATCAGCATA
TGACCCCCGTTTA
>Zotu10
CGTCGAAAAATTT
CGATACACCCTAT
>Zotu22
CGTACGTCCCCTT
CGATATAATATATA

$ cat test.txt
Zotu1
Zotu2
Question: Hi all,

I have Illumina sequencing results of a bacterial genome and a quality score of 35.89 is associated with these data. I know that a quality score of 30 is 99.99% of base calling accuracy based on [this][1] but what about the meaning of 35.89? Is there a scale used for quality score?

Thank you very much for your much appreciated help on this!

[1]: https://www.illumina.com/documents/products/technotes/technote_Q-Scores.pdf

Answer: The phred score describes the probability of a base call being incorrect. [Wikipedia](https://en.wikipedia.org/wiki/Phred_quality_score) has a great explanation including the formula, which is `-10xlog10(probability of being wrong)`.

```
phred = 35.89 # = -10*log10(probability)
35.89/-10 = -3.589
probability = 10^-3.589 = 0.0002576321
```
Therefore, a base call with a Phred score of 35.89 has a 0.02% chance of being wrong. Or, in other words, 0.02% of bases with a Phred score of 35.89 will have been wrongly assigned.
Question: Hi,

I'm trying to use Batch Entrez (https://www.ncbi.nlm.nih.gov/sites/batchentrez) to look up the identifiers below, but although UID are recognised (and return records in the CBI database) - the output is empty.

Any idea why?
Thanks! Reproducible example (to copy paste in a .txt document):

KF771749.1
KF153813.1
MG737737.1
MK091921.1
MG173246.1
JQ073586.1
HM402114.1
MG851790.1
KP317997.1
MG737738.1

Answer: Using EntrezDirect:

$ esearch -db nuccore -query "MG737737" | elink -target taxonomy | esummary | xtract -pattern DocumentSummary -element TaxId | xargs -n 1 sh -c 'efetch -db taxonomy -id "$0" -format xml' | xtract -pattern Taxon -block "*/Taxon" -unless Rank -equals "no rank" -tab "\n" -element Rank,ScientificName
superkingdom	Eukaryota
clade	Opisthokonta
kingdom	Metazoa
clade	Eumetazoa
clade	Bilateria
clade	Protostomia
clade	Ecdysozoa
clade	Panarthropoda
phylum	Arthropoda
clade	Mandibulata
clade	Pancrustacea
subphylum	Hexapoda
class	Insecta
clade	Dicondylia
subclass	Pterygota
infraclass	Neoptera
cohort	Paraneoptera
order	Hemiptera
suborder	Auchenorrhyncha
infraorder	Cicadomorpha
superfamily	Cicadoidea
family	Cicadidae
subfamily	Cicadettinae
tribe	Cicadettini
genus	Kikihia
Question: I want to delete the part starting after ">" until "cds_". And the characters after the accession. In the example "_1" and "_2".
In the original file this counts up to 784, so "_784". Can someone help me with a solution? Would be great.

From this

>lcl|NC_002712.2_cds_NP_268424.1_1 [gene=dnaZA]
ATGACTG
>lcl|NC_002721.2_cds_NP_268453.1_2 [gene=dnaGC]
ATGTTCG

To

>NP_268424.1 [gene=dnaZA]
ATGACTG
>NP_268453.1 [gene=dnaGC]
ATGTTCG

Answer: cat $filename | sed 's/>.*cds_/>/' | sed -e 's/_[0-9]* \[/ \[/'

*First sed*

Each character (.*) that is contained between `>` and `cds_` is replaced by only `>`. This leaves the FASTA header intact but removes everything until `NP_`.

*Second sed*

Enabling regular expressions (`-e`) we substitute each combination of underscore (`_`), a series of numeric characters (`[0-9]*`), a whitespace (` `) and a bracket, escaped (` \[`) by a whitespace and an escaped bracket only (` \[`).
Question: I've noticed that CellRanger defaults to running STAR on sc-RNA-seq samples considering only reads that align to the exome.

However, they suggest modifying the default when running nuclei samples to align to exons + introns (pre-RNA).

Given that single-cell data includes cytosolic + nuclear RNA, why is it not standard practice to use the pre-RNA alignment strategy as the default for sc-RNA-seq experiments?

Since it generally improves alignment, it seems like a good idea, but am guessing there may be unintended consequences I'm not aware of.

Thanks in advance for any insights you might have!

Answer: I received this response after inquiring to 10x directly (summarized):

Other mechanisms could also contribute to intronic reads, such as internal poly-A priming (see the detailed [technical note][1]). By including intronic read counting, it will potentially have a bias for poly-A-rich genes (not referring to exons here). Aligning with introns increases the sensitivity of gene detection, since it helps the detection of genes that have intronic UMIs but 0 exonic UMIs in any cell.

Different samples have different intronic read mapping profiles (table 1 in above link). (Also, if there is a polyA stretch within an exon, then multiple 'unique' reads could potentially be generated from a single transcript). Therefore, if you expect a high fraction of intronic reads in your single-cell sample and would like to increase the gene detection sensitivity, you could include --include-introns.

[1]: https://assets.ctfassets.net/an68im79xiti/awNZTarmwqmwxKcvDv9wv/b05500661e36290b8ce59689ff889ea8/CG000376_TechNote_Antisense_Intronic_Reads_SingleCellGeneExpression_RevA.pdf
Question: library(GenomicRanges)

gr <- GRanges(Rle(c("chr2", "chr2", "chr1", "chr3"), c(1, 3, 2, 4)),
IRanges(1:10, width=10:1, names=head(letters, 10)),
Rle(strand(c("-", "+", "*", "+", "-")), c(1, 2, 2, 3, 2)),
score=1:10, GC=seq(1, 0, length=10))
split(gr, seqnames(gr))

Gives me a list of chromosomes, but how do I get a list split on chromosome and strand?

Answer: Paste the chrom and strand then split as usual:

split(gr, paste(seqnames(gr), strand(gr)))