cannin/biostars_qa huggingface dataset:
Question: I am advised that the FASTA sequences in Uniprot is the wild type. However when I check the sequence of HIV-1 protease in Uniprot, I see it isn't the same as in 1MUI, the wild type as this side stated.

The FASTA in Uniprot:

PQVTLWQRPIVTIKIGGQLKEALLDTGADDTVLEEMSLPGKWKPKMIGGIGGFIKVRQYDQVSIEICGHKAIGTVLIGPTPVNIIGRNLLTQLGCTLNF

The one in 1MUI:

PQITLWQRPLVTIKIGGQLKEALLDTGADDTVLEEMSLPGRWKPKMIGGIGGFIKVRQYDQILIEICGHKAIGTVLVGPTPVNIIGRNLLTQIGCTLNF

Which is the true wild type?

Answer: I am not a HIV expert, but as far as I know there is so much variability in the HIV-1 protease, that there is no such thing as a wildtype or consensus sequence for this protein.

Also, it is good to realise that UniProtKB consists of two parts: UniProtKB/Swiss-Prot, which is manually annotated and reviewed and therefore of very high quality, and UniProtKB/TrEMBL, which is automatically annotated and not reviewed and therefore of much lower quality. See also this FAQ: ["Why is UniProtKB composed of 2 sections, UniProtKB/Swiss-Prot and UniProtKB/TrEMBL?"][1].

[1]: http://www.uniprot.org/faq/7
Question: Hi all,

I want to use enhancedvolcanoes package to produce volcanoes plot, but the example is so confusing with its input data in 'S4' type, obj oriented. I have been looking for an example that shows simple data format that I need to have before inputing them in and use the package. Nowhere to be found!

I just have an excel sheet with 3 columns: `Gene, log2FoldChange, pvalue`

When input into R as VC_data, can we seen as:

https://www.dropbox.com/s/pq1l38ko1p35e8k/Screen%20Shot%202019-05-19%20at%2011.22.13%20PM.png?dl=0

My code is as following:

EnhancedVolcano(VC_data,
lab = VC_data$Gene,
x = 'log2FoldChange',
y = 'pvalue',
xlim = c(-5, 8))

I keep getting this error saying that

'Error in EnhancedVolcano(VC_data, lab = VC_data$Gene, x = "log2FoldChange", :
log2FoldChange is not numeric!'

When I did `is.numeric(VC_data$log2FoldChange)`, I get true.

I am very frustrated. Please help!

Answer: Seems to be an issue with your input being a `tibble`:

tb <- tibble(log2fc = c(2,3), pval = c(0.01, 0.007))

## this returns the error you have:
EnhancedVolcano(toptable = tb, x = 'log2fc', y = 'pval', lab = "foo")

## this one does it properly:
EnhancedVolcano(toptable = data.frame(tb), x = 'log2fc', y = 'pval', lab = "foo")

=> Convert `tibble` to `data.frame`

Edit:
The way EnhancedVolcano checks if data are numeric is (source code line 76-78):

if(!is.numeric(toptable[,x])) {
stop(paste(x, " is not numeric!", sep=""))
}

Problem is that e.g. `tb[,1]` does not automatically reduces dimensions to directly access the data but returns:

> tb[,1]
# A tibble: 2 x 1
log2fc

1 2
2 3

whereas:

tb[,1][[1]]
[1] 2 3

returns the numeric data.
Question: Hi guys,

I am looking for a solid tool for viral metagenomic data analysis, for virus profiling etc. I am new in this domaine, I read some papers which propose several tools, I need to pick one and stick on it. So It is better ask the community first. I found tool such as **Virsorter**, **VirFinder**,**PPR-Meta**,**DeepVirFinder** and **VirionLang**, have no idea which one is relatively the best.

Please help.

My Metagenomics data is from human fecal sample, medium depth. I've already run humann3 for shallow microbiome analysis. Yes, humann3(or metaphlan3) can do virus profiling with additional parameter `--add_viruses`, but results was not good enough.

Answer: I don't think there is a best tool for this purpose, as all of them them are better at one thing or another. Generally speaking, homology-based tools tend to have lower false positive rate, yet tools that rely on k-mer frequencies may work better in discovering novel and rare viruses. The most reliable solution is likely to come as a union of several approaches.

I recommend [**VIBRANT**][1] as a standalone tool. Yet another approach is a combination of VirSorter2, CheckV and DRAM-v as detailed in [**this protocol**][2] by Sullivan lab.

I don't think you need pleas for help in your messages. It is understood that everyone who writes in this forum needs help.

[1]: https://github.com/AnantharamanLab/VIBRANT
[2]: https://www.protocols.io/view/viral-sequence-identification-sop-with-virsorter2-bwm5pc86
Question: Common Workflow Language (**CWL**) https://github.com/common-workflow-language/common-workflow-language / http://common-workflow-language.github.io/draft-3/ has been trending on my twitter timeline during the last weeks.

However the spec is quite large and I find it hard to get some simple examples.

Furthermore, I have the feeling that all engines require a lot of dependencies or docker. I'd like to test my makefile-based workflows using CWL, how should I write and test the following simple **Makefile** using CWL:

```
SHELL=/bin/bash
.PHONY: all
all : database.dna

database.dna : seq1.dna seq2.dna seq3.dna
cat seq1.dna seq2.dna seq3.dna > database.dna

seq3.dna : seq3.rna
tr "U" "T" < seq3.rna > seq3.dna

seq3.rna :
echo "AUGCGAUCGAUCG" > seq3.rna

seq2.dna : seq2.rna
tr "U" "T" < seq2.rna > seq2.dna

seq2.rna :
echo "AUGAAGACUGCGAUCGAUCG" > seq2.rna

seq1.dna : seq1.rna
tr "U" "T" < seq1.rna > seq1.dna

seq1.rna :
echo "AUGAAGACUGACUCGUCG" > seq1.rna
```

**EDIT**: feel free to add the file for your favorite workflow-engine as an answer.

https://twitter.com/PaoloDiTommaso/status/625995681434607616

https://twitter.com/smllmp/status/625999447869231104

Answer: To describe this workflow using CWL you need to:

- describe the 3 different tools (tr, echo, cat)

- describe the workflow structure

- describe the input data for the workflow job, i.e. the 3 input strings (AUGAAGACUGACUCGAUCGAUCG. etc.)

As an example to see the implementation of "tr", you can see this gist: https://gist.github.com/hmenager/897f3e81ad7cd98e94a6 . The github repository for common-workflow-language also has a number of tool and workflow descriptions that might help (https://github.com/common-workflow-language/common-workflow-language/tree/master/conformance/draft-2).

The reference implementation should be fairly easy to set up (pip install cwl-runner) but it is minimal (no cluster integration, etc.). Docker is not and should never become a requirement of CWL, some tools have been described with docker "requirements" but this is absolutely not mandatory. You can also ask your questions on all the channels mentionned here: https://github.com/common-workflow-language/common-workflow-language#community-and-contributing

Hope this helps,

Hervé
Question: I'm working on a SLURM cluster with NGS data. I trimmed raw reads and was thinking of the best way to align them to the reference genome. I have pairs of reads for a few samples. I wrote a script for parallel bwa:

#SBATCH --cpus-per-task=1
#SBATCH --ntasks=10
#SBATCH --nodes=1

# align with bwa & convert to bam
bwatosam() {
id=$1
index=$2
output=$3/"$id".bam
fq1=$4/"$id".R1.fq.gz
fq2=$4/"$id".R2.fq.gz

bwa mem -t 16 -R '@RG\tID:"$id"\tSM:"$id"\tPL:ILLUMINA\tLB:"$id"_exome' -v 3 -M $index $fq1 $fq2 |
samtools view -bo $output
};
export -f bwatosam

# run bwatosam in parallel
ls trimmed/*.R1.fq.gz |
xargs -n 1 basename |
awk -F ".R1" '{print $1 | "sort -u"}' |
parallel -j $SLURM_NTASKS "bwatosam {} index.fa alns trimmed"

But I'm not sure if I use the right parameters (#SBATCH) for the job because if I do it without -j:

#SBATCH --nodes=1
#SBATCH --ntasks-per-node=5

# run bwatosam in parallel
ls trimmed/*.R1.fq.gz |
xargs -n 1 basename |
awk -F ".R1" '{print $1 | "sort -u"}' |
parallel "bwatosam {} index.fa alns trimmed"

It works 10 times faster. What number of nodes/cpus/threads should I use?

Answer: Depends on the node. I typically run alignments with basically this kind of script (sorry https://www.biostars.org/u/30/, no snakemake yet) on a 72-core node with 192GB RAM, and then use:

#SBATCH --nodes=1
#SBATCH --ntasks-per-node=72
#SBATCH --partition=normal

In this case I would use 4 parallel processes with 16 threads for `bwa` each. Depends on how much memory your node has. Can you give some details? When using `parallel` I recommend booking the entire node to ensure you are not interfering with processes from other users.

=> Note that I always book the entire node if running `parallel` things so I essentially do not care about RAM consumption etc. as long as the node can handle it. If you share the node with others it might be a good idea to task the your admin before if using `parallel` is allowed on your cluster nodes.
Question: I am trying to understand Gene expression data analysis. Could someone the percentage of change of the GE calculating from log2 fold change (L2FC)?

> FoldChange = μ_Treatment / μ_Control
>
> log2FC=Log2(μ_Treatment)-Log2(μ_Control)

Answer: seems correct to me. log2FC is more informative as you will see the direction (sign) and the effect.

Treatment = 10
Control = 100

`FoldChange = 10/100 = 0.1` which is not informative

`log2FC = log2(10)-log2(100) = -3.3` which tells the effect and direction.
Question: I have a bed file with the co-ordinates includes chr start and end columns.

#chr	start	end
chrM	1	1
chrM	2	2
chrM	3	3
chrM	4	4
chrM	5	5
chrM	7	7
chrM	8	8
chrM	9	9
chrM	10	10
chrM	11	11
chrM	12	12
chrM	13	13
chrM	14	14
here there are 2 region are present i want to get output like

#chr start end
chrM 1 5
chrM 7 14
which tool should i use to perform this
after this i need to find its location into the mitochorial genome

Answer: I'm afraid your BED file only describes **deletions** because a BED in half-open : see chromEnd in https://genome.ucsc.edu/FAQ/FAQformat.html#format1

the tool you need is "bedtools merge" http://bedtools.readthedocs.io/en/latest/content/tools/merge.html

> which tool should i use to perform this after this i need to find its location into the mitochorial genome

uh ???
Question: Hello,

I merged bam files using `samtools merge -r out.bam in1.bam in2.bam ...`

With `-r` option I got the RG tag for each read, as expected, good. However, the header of the merged bam does not have the @RG lines. So my question: Is there any off-the-shelf tools to add the @RG header lines once the reads have been tagged?

If not, the pipeline I have in mind goes along these lines:

- Scan the bam file and collect all the different RG tags
- Output the header of the bam file.
- Append the RG tags to the header the
- Reheader the original bam file.

Does it make sense?

Thanks

Dario

Answer: I ended up writing a little script for this, if anyone is interested it is here:

- https://code.google.com/p/bioinformatics-misc/source/browse/trunk/addRGtoSAMHeader.py
- https://github.com/dariober/bioinformatics-cafe/blob/master/addRGtoSAMHeader.py
Question: Hi everyone,

I came across SNP filtering tutorial where the author used the flag `--mac 3` to filter SNPs that have a minor allele count less than 3. That is `vcftools --vcf input_file.vcf --mac 3 --recode --out filtered_file`
Could someone explain to me why filtering out sites with minor allele count below 3? By retaining 3 alleles and above, what exactly are we aiming at?
I tried to apply the above script to my snps data from cassava crops (diploid with 18 chromosomes) having 359793 sites x 980 samples. After filtering, I now have 147518 sites x 980 samples indicating large number of sites where dropped.
What is going on? Please I need better clarification on --mac 3. This is because I also intend to filter for minor allele frequency later on.enter code here
Thanks

Answer: Alleles which are observed infrequently are more likely to be errors and in any case, they would not be useful because they have low statistical power to detect association with anything so they are usually filtered out.
Question: H!
Which bioinformatics tools supporting CRISPR-Cas Experiments do you recommend? Or which tools are you using?

Answer: Here is [CHOPCHOP][1] which has been developed at our institute. It supports different approaches (Cas9, Cas13, Cpf1, TALEN) and knock-out, knock-in, activation, and repression.

And here is a more comprehensive list of all [CRISPR tools in bio.tools][2] sorted by citations, as in "most used successfully".

[1]: https://chopchop.cbu.uib.no/
[2]: https://bio.tools/t?q=crispr&sort=citationCount&ord=desc
Question: Hi All,

Is it possible to re-run STAR RNA aligner to obtain wiggle plots on already aligned and sorted BAM files? Or do I need to re-align again to obtain this type of output? I really do not want to have to re-align - so I am open to other options.

Thank you!

Answer: See https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4631051/

STAR --runMode inputAlignmentsFromBAM \
--inputBAMfile Aligned.sortedByCoord.out.bam \
--outWigType bedGraph \
--outWigStrand Stranded
Question: Hi,
I would like to calculate the percent-identity from a CIGAR string from a BAM/SAM file containing alignments. I want to calculate the PID only for the aligned region, ignoring clipped ends ("H","S"). I can parse the CIGAR in R and get the sums for each letter in the CIGAR, so it's just a conceptual question about the definition and whether or not it is correct:
so given the CIGAR contains characters: "M" (match),"N" (skip),"D" (Deletion),"I"(Insertion), "S" (soft-clip),"H" (hard-clip), of which I ignore "S", "H", and "M" is the total sum of M, N the total of N, etc.:

e.g.: 10S 20M5I5D20M 10S

Is this a good way of defining the formula?

pid1 := 100 * M / (M+N+I)

or maybe

pid2 := 100* M / (M+N+I-D)

the example: pid1 = ~88% but then pid2 = 100% (which wouldn't make much sense).

Related: http://biostar.stackexchange.com/questions/9358/is-there-any-r-package-to-parse-cigar-element-of-sam-files/17031#17031

Thank you very much

Answer: I would say pid1 would be correct. Computationally, I don't think there is any difference between N (skipped) and D (deletion).

Both skipped and deletion should produce something like this:

read ACGTACGT--ACGTACGT
reference ACGTACGTAAACGTACGT

So if you add N + D, you are adding the gap twice.

edit **

The SAM specs is a bit confusing on the matter:

M alignment match (can be a sequence match or mismatch)
I insertion to the reference
D deletion from the reference
N skipped region from the reference
S soft clipping (clipped sequences present in SEQ)
H hard clipping (clipped sequences NOT present in SEQ)
P padding (silent deletion from padded reference)
= sequence match
X sequence mismatch

M can be a match or mismatch??

Further down in the sam specs:

For mRNA-to-genome alignment, an N operation represents an intron. For other types of alignments, the interpretation of N is not defined.

Sum of lengths of the M/I/S/=/X operations shall equal the length of SEQ
Question: I have gtf file containing 7 columns from which i am going extract only gene name part from column 7th.

The 7th column contains some information including `gene_id`, `gene_name` and so on that I posted one row of column 7th below:

gene_id "XLOC_000001"; transcript_id "TCONS_00000001"; exon_number "1"; gene_name "NAC001"; oId "AT1G01010.1"; nearest_ref "AT1G01010.1"; class_code "="; tss_id

I need only `gene_name` from this column, for example "NAC001", how can I extract which I need please?

Answer: Assuming the fields in your gtf are tab separated you could try something like this:

awk 'BEGIN{FS="\t"}{print $7}' YOURFILE | awk 'BEGIN{FS="gene_name"}{print $2}' | awk 'BEGIN{FS=";"}{print $1}' > OUTFILE

genomax2 solution with sub is nicer.
Question: I made a script to retrieve all of the recorded full-length human lung infecting coronaviruses in genbank (virus variation database, ncbi). After trimming out all of the really small files or those that didnt work, I'm left with ~200 fasta files to work with. What are some analyses I can run with the new COVID-19 genome? I made a blast database and ran some blastn queries but I'm wondering what bioinformatic analyses are typically run on novel viruses?

My blastn query of word size 7 had a few hits but I am wondering how to interpret these results and what to do with that data.

This is just for fun/educational. I don't have much experience with viruses or comparative genomics.

Thanks

Answer: Since you selected these strains (which are similar) there is not much of a point in doing BLAST analysis. Instead, you can start doing a multiple sequence alignment with the sequences. If you are serious about learning then there are command line version of `MAFTT`, `T-COFFEE`, `Clustal` , `MUSCLE`. You will also find web front-ends for these tools (if you search around). A multiple sequence alignment gives you an idea of relationship of the sequences to each other. This can then be used to infer possible evolutionary relationships among the genomes.

`MEGA` is [user friendly software][1] that has a GUI for doing above analysis. They have a [pretty good online][2] manual as well.

[1]: https://www.megasoftware.net/
[2]: https://www.megasoftware.net/web_help_10/index.htm#t=Preface.htm
Question: Hi, I have two multifasta files. I want to merge them deleting all those fasta seqences from the first multifasta file which are also in the second file. I need to do it by header comparison, sequences are different under the same headers.

Alternatively, could somebody give me a hint how to generate all the contigs (even those unchanged) through bcftools consensus?

Thanks,
Pawel

Answer: You can use the BBMap package like this:

filterbyname.sh in=file1.fasta names=file2.fasta exclude out=file1_filtered.fasta
cat file1_filtered.fasta file2.fasta > combined.fasta
Question: I am trying to assembled a pooled metagenomic dataset with 762909004 reads. Will using tools such as BBNorm to downsize the coverage in order to speed up the assembly significantly effect the final assemblies?

Answer: To answer your question first: normalizing it down to 100x should have no ill effect, and it may even help reduce the errors if you have very high coverage. Several times I have obtained better assemblies from 60x- or 80x-normalized data than from original datasets. Still, I suggest that you assemble from your original data so you have a baseline. I have assembled about 350 million reads with [**megahit**][1] in about 10 hours - you should be able to do it with your dataset in couple of days. It is always possible the original data is of such quality that it yields the best assembly.

[1]: https://github.com/voutcn/megahit
Question: Good Morning,

I'm a student of Bioinformatics from the University of Bologna.

I'm working on a personal idea and I'm using MUSTANG to perform multiple structural alignment.
To run the algorithm .pdb file with just one chain are required so I'm looking for tools or a suggestion in order to edit in a clever and fast way all the files.

'til now I wrote some command lines in python but I think that they're not enough specific to correctly edit all the files in fact I often have to check manually all the structures again.

someone could help me?

thanks in advance!
Dade.

Answer: """
Extract a single chain from a PDB file
"""

from __future__ import print_function
import Bio.PDB
import Bio.PDB.PDBIO
import sys
import argparse

class ChainSelect(Bio.PDB.Select):
def __init__(self, target_chain):
self.target_chain = target_chain
def accept_chain(self, chain):
if chain.get_id() == self.target_chain:
return 1
else:
return 0

def main():
argparser = argparse.ArgumentParser(description="Extract chain from a PDB file")
argparser.add_argument('infile', help="Path to input file (PDB)")
argparser.add_argument('chain', help="Chain to extract")
argparser.add_argument('outfile', help="Path to output file (PDB)")
args = argparser.parse_args()

pdbparser = Bio.PDB.PDBParser()
io = Bio.PDB.PDBIO()
with open(args.infile, 'r') as infile:
struct = pdbparser.get_structure(args.infile, infile)
io.set_structure(struct)
with open(args.outfile, 'w') as outfile:
io.save(outfile, ChainSelect(args.chain))
return 0

if __name__=="__main__":
sys.exit(main())

Run it with something like:

python extract_chain.py 1XXX.pdb A 1XXXA.pdb

to extract just chain A from 1XXX.pdb.
Question: Hi all

This question could be considered as a follow up of this discussion.
http://www.biostars.org/post/show/3407/how-to-extract-reads-from-bam-that-overlap-with-specific-regions/#3414
What I need is to extract reads from bam file that fall only within a given region (not overlap the given region), the region being in the form of a gff file or bed file. Overlapping reads could be extracted by several methods (as in the discussion mentioned or BEDTools). The idea is to try to be pretty sure of excluding 5´ UTRs in the process of detecting intergenic transcripts. I saw a tool in BamUtil (http://genome.sph.umich.edu/wiki/BamUtil) called "writeRegion" which would pretty much do what I want. Somehow could not get this running for my dataset.
Was wondering if you guys might have an "R" or some other solution for this.
Thanks in advance
Abi

Answer: You can extract mappings of a sam/bam file by reference and region with samtools. For example:

samtools view input.bam "Chr10:18000-45500" > output.bam

That would output all reads in Chr10 between 18000-45500 bp.
Question: I am trying to modify the location of features within a GenBank file. I know `feature.type` will give gene/CDS and `feature.qualifiers` will then give "db_xref"/"locus_tag"/"inference" etc. Is there a `feature.` object which will allow me to access the location (eg: `[5240:7267](+)`) directly?

This URL give a bit more info, though I can't figure out how to use it for my purpose... http://biopython.org/DIST/docs/api/Bio.SeqFeature.SeqFeature-class.html#location_operator

Essentially, I want to modify the following bit of a GenBank file:

gene 5240..7267
/db_xref="GeneID:887081"
/locus_tag="Rv0005"
/gene="gyrB"
CDS 5240..7267
/locus_tag="Rv0005"
/inference="protein motif:PROSITE:PS00177"
...........................

to

gene 5357..7267
/db_xref="GeneID:887081"
/locus_tag="Rv0005"
/gene="gyrB"
CDS 5357..7267
/locus_tag="Rv0005"
/inference="protein motif:PROSITE:PS00177"
.............................

Note the changes from **5240** to **5357**

So far I have the following python script:

from Bio import SeqIO
gb_file = "mtbtomod.gb"
gb_record = SeqIO.parse(open(gb_file, "r+"), "genbank")
rvnumber = 'Rv0005'
newstart = 5357

final_features = []

for record in gb_record:
for feature in record.features:
if feature.type == "gene":
if feature.qualifiers["locus_tag"][0] == rvnumber:
if feature.location.strand == 1:
# Amend feature location from current to 'newstart'
else:
# do the reverse for the complementary strand
final_features.append(feature)
record.features = final_features
with open("test.gb","w") as test:
SeqIO.write(record, test, "genbank")

*Rv0005* is just an example of a locus_tag I need to update. I have about 600 more locations to update per GenBank file, and about 10-20 GenBank files to process (with more to come)

Answer: Where is the question? Do you get an error message?

Well I do get one: when I do something like:

features.location.start.position = newstart

I get an error:

seqfeature AttributeError can't set attribute

And it doesn't matter whether the genbank file was open read-only as you did, or in writing mode.

I suppose that you need to create a new genbank file, and copy all the attributes.

**Edit 1**:

Well modifying the start is possible after all.

```
from Bio import SeqIO
from Bio import SeqFeature
start_pos = SeqFeature.AfterPosition(newstart)
end_pos = SeqFeature.BeforePosition(feature.location.end.position)
my_location = SeqFeature.FeatureLocation(start_pos, end_pos)
features = my_location
```

But that does not solve the saving problem yet.

**Edit 2**:

Useful links:

- https://www.biostars.org/p/57549/
- [Dealing with GenBank files in Biopython][1]
- [Parsing Genbank files with Biopython][2]

[1]: http://www2.warwick.ac.uk/fac/sci/moac/people/students/peter_cock/python/genbank/
[2]: http://wilke.openwetware.org/Parsing_Genbank_files_with_Biopython.html