cannin/biostars_qa huggingface dataset:
Question: Hi everyone,

I am trying to found out what residues are in contact with a ligand (the binding site) in a PDB file. To give you an example, I need to find the residues in contact with FMN ligand in the PDB entry 1AL7.

First of all, can we find this type of information in a PDB file? Or do we need a software to calculate it?
I think the first question is yes, but I am not sure: I used the new mmCIF format, as they are well organized compared to the old version. In this format, there is a Data Category called "**chem_comp**" and the dictionnary on the PDB website says :

"*Data items in the CHEM_COMP category give details about each of the chemical components from which the relevant chemical structures can be constructed, such as name, mass or charge.*"
It allows me to know if there is the presence of a FMN ligand.

Another category is "**struct_site**":

*"Data items in the STRUCT_SITE category record details about portions of the structure that contribute to structurally relevant sites (e.g. active sites, substrate-binding subsites, metal-coordination sites)."*

In my example, there is a line "AC1 Software ? ? ? ? 22 'BINDING SITE FOR RESIDUE FMN A 360'"

And then, you have another category called "**struct_site_gen**":

"*Data items in the STRUCT_SITE_GEN category record details about the generation of portions of the structure that contribute to structurally relevant sites.*"

This category is, I guess, my answer, because I see all the residues with their number sharing the same id (AC1) than in struct_site.

Can I use this information for the binding site or you guys think I have to use a software to calculate the residues in the binding site? If so, which one?

Thank you

Answer: Yes, you can use the information in STRUCT_SITE_GEN listed under AC1 to get the binding residues: TYR24, TYR25, ALA76 etc.

If you did want to do this from the coordinate data, one way would be to use [BioStructures.jl][1] in Julia:

using BioStructures
struc = read("1AL7.cif", MMCIF)
res = collectresidues(struc, standardselector)
for r in res
if distance(r, struc["A"]["H_360"]) < 4.0
println(r)
end
end

This returns more or less the same residues as STRUCT_SITE_GEN. You could also do this with Biopython.

[1]: https://github.com/BioJulia/BioStructures.jl
Question: Hi I have a question regarding when steps are able to run in parallel. We noticed that sometimes all of the jobs from a scatter must complete before moving to the next

I was wondering whether a scattered subworkflow can start before all the inputs to that step are available, or does a subworkflow have to wait for all of the previous inputs to be available regardless of whether they are available for **just the first set** of scattered inputs?

Answer: All inputs for the step have to be available. See [this][1], look for bullet point "Instead of scattering separate steps, prefer to scatter over a subworkflow."

[1]: https://doc.arvados.org/user/cwl/cwl-style.html
Question: Any tool that can do median per base coverage over pre-defined windows from a BAM file?

I could make the actual per base coverage using `samtools depth` or `genomeCoverageBed`.

chr pos read_depth
1 1 6
1 2 8
...

I already have the windows prepared using

bedtools makewindows \
-g "genome.fa.fai" \
-w 1000 > "windows.bed"

I am looking for a result like below:

Median per base coverage over 1000bp windows

chr start end median_read_depth
1 1 1000 5
1 1000 2000 8
...

Answer: I 've written : http://lindenb.github.io/jvarkit/BamStats04.html

$ java -jar dist/bamstats04.jar -B src/test/resources/toy.bed.gz src/test/resources/toy.bam 2> /dev/null | column -t

#chrom start end length sample mincov maxcov meancov mediancov nocoveragebp percentcovered
ref 10 13 3 S1 3 3 3.0 3.0 0 100
ref2 1 2 1 S1 2 2 2.0 2.0 0 100
ref2 13 14 1 S1 6 6 6.0 6.0 0 100
ref2 16 17 1 S1 6 6 6.0 6.0 0 100
Question: Is there an option in Hisat that is similar to `-g`/`--max-multihits` in Tophat? I can not find it in the parameter list on https://ccb.jhu.edu/software/hisat/manual.shtml#options. There is a Uniqueness class in here: https://github.com/infphilo/hisat/blob/master/unique.h, but when I grep for that class or the `bestIsUnique` function in the Hisat source folder I can't find it used anywhere.

Or would it work to use `samtools view -q 1`?

If anyone else is interested, I e-mailed the author and he said:

> One possible solution is to use the next best alignment score (ZS:i:, previously XS:i:). If the next alignment score is the same as the primary alignment score (AS:i:), that means there are at least two equivalent alignments found by HISAT. This additional SAM field (ZS:i:) is only available at HISAT github, which I'll include in the next release of HISAT.

Answer: The most recent release of hisat (v0.1.6) uses NH tags in the BAM to specify the number of alignments for the read/pair (see release notes: https://ccb.jhu.edu/software/hisat/index.shtml). As far as I understand it, reads/pairs with NH:i:1 have a single best alignment and are therefore "uniquely" mapped.
Question: I'm using ete toolkit to get a phylogenetic tree for a list of ncbi taxids:

```py
from ete2 import NCBITaxa
ncbi = NCBITaxa()
tree = ncbi.get_topology([9606, 9598, 10090, 7707, 8782])
print tree.get_ascii(attributes=["sci_name", "rank"])
```

Printing it with ascii chars works, but how to render this tree *including* the attributes (sci_name, rank) to an image?

tree.render("tree.pdf")

There seems to be no 'attribute' in this function.

Answer: `tree.render()` is a general purpose method. As node attributes are completely arbitrary, you need to specify what should be drawn and where... Check the docs regarding the drawing system: https://pythonhosted.org/ete2/tutorial/tutorial_drawing.html

For your example, something like this should work:

```py
from ete2 import NCBITaxa, AttrFace, TreeStyle
ncbi = NCBITaxa()
tree = ncbi.get_topology([9606, 9598, 10090, 7707, 8782])

# custom layout: adds "rank" on top of branches, and sci_name as tip names
def my_layout(node):
if getattr(node, "rank", None):
rank_face = AttrFace("rank", fsize=7, fgcolor="indianred")
node.add_face(rank_face, column=0, position="branch-top")
if node.is_leaf():
sciname_face = AttrFace("sci_name", fsize=9, fgcolor="steelblue")
node.add_face(sciname_face, column=0, position="branch-right")
ts = TreeStyle()
ts.layout_fn = my_layout
ts.show_leaf_name = False
tree.render("tree.png", tree_style=ts)
```
Question: Hi folks!

I have two GFF3 files, one with annotated genes. Another with predicted transposons.

I want to find the numerical distribution of the distances between transposons and genes.

I already found the distance amongst just transposons, and amongst just genes - using the 'spacing' sub-command of bedtools 2.2.4.0, however 'spacing' cannot compare two GFF files.

For comparing spacing of features in 2 separate GFFs, bedtools 2.24.0 has the reldist sub-command, but this does **not** yield absolute distance in base pairs, but only relative distance distribution, so its not so useful to me...

So I am not sure if there is an off the shelf option in bedtools that can help me answer this question. Would any of you have a simple solution to my problem? It does not have to be using bedtools... Thank you!

Answer: Maybe `closestBed` with `-d` option is what you want?

```
Tool: bedtools closest (aka closestBed)
Version: v2.23.0
Summary: For each feature in A, finds the closest
feature (upstream or downstream) in B.
...
-d In addition to the closest feature in B,
report its distance to A as an extra column.
- The reported distance for overlapping features will be 0.
```
Question: What happened with miRecords mirna database website? The webpage is not available and I used it's data as part of an analysis in my paper. Now the reviewer can't see the page and asks for an explanation.

Edit:

Sorry, I forgot to post the link. The miRecords was suppose to be available at http://miRecords.umn.edu/miRecords.

Answer: http://mirecords.umn.edu/miRecords/ works fine for me at this moment. It may just have been a hiccup.
Question: Hi all,

I'm dealing with a fasta file with spaces at the end of line, which caused the problem. I didn't find a suitable way to remove them. Please kindly tell me the appropriate command for removing them?

Answer: Oh my gawk!

All previous solutions would risk modifying your fasta header as well. This one will not.

gawk 'BEGIN{line=0}{ if ($0 !~/^>/ && $0 ~/ +/ ) {gsub(/ +/, //); line++} print}END{print line" lines with white spaces treated" > "/dev/stderr"}' myfasta.fa >output.fa

If you only want to remove the spaces at the end of the lines:

gawk 'BEGIN{line=0}{ if ($0 !~/^>/ && $0 ~/ +$/ ) {gsub(/ +$/, //); line++} print}END{print line" lines with white spaces treated" > "/dev/stderr"}' myfasta.fa>output.fa
Question: Hello everyone (first post ever for me)

I am using Biomart in R and I am trying to obtain the maize homologs of a set of barley genes (I am only displaying one here for example's sake). Everything works fine when I use the getBM function for getLDS does not seem to work.

```r
host="plants.ensembl.org"
mysets<-listDatasets(useMart("plants_mart", host = host))
barley <- useDataset("hvulgare_eg_gene", mart = useMart("plants_mart", host = host))
zea <-useDataset("zmays_eg_gene", mart = useMart("plants_mart", host = host))
getLDS(attributes = c("ensembl_gene_id","description", "go_id","name_1006"),
filters = "ensembl_gene_id", values = "HORVU.MOREX.r3.1HG0000220",
mart = barley ,
attributesL = c("ensembl_gene_id","description", "go_id","name_1006"),
martL = zea)
```

I get the following error message. I think the issue may be that the Virtual Schema is set to default, but I can't find how to change it to plants_mart

```r
Error in getLDS(attributes = c("ensembl_gene_id", "description", "go_id", : Query ERROR: caught BioMart::Exception::Usage: WITHIN Virtual Schema : default, Dataset hvulgare_eg_gene NOT FOUND
Traceback:

1. getLDS(attributes = c("ensembl_gene_id", "description", "go_id",
. "name_1006"), filters = "ensembl_gene_id", values = "HORVU.MOREX.r3.1HG0000220",
. mart = barley, attributesL = c("ensembl_gene_id", "description",
. "go_id", "name_1006"), martL = zea)
2. stop(postRes)
```

Any help would be greatly appreciated

Answer: Great spot that this was related to the `virtualSchemaName`. That was hardcoded to be `"default"` in **biomaRt** and no one has ever reported a problem before. I assumed there were no Mart instances that ever used anything different!

Anyway, it's been patched in **biomaRt** version 2.50.3, which should be available from Bioconductor in the next few days. If you want to get a copy before then you can install from Github with `BiocManager::install('grimbough/biomaRt', ref = 'RELEASE_3_14')`.

Here's an example with it working for a similar query:

```r
library(biomaRt)
packageVersion("biomaRt")
#> [1] '2.50.3'

ensembl_plants <- useEnsemblGenomes(biomart = "plants_mart")

barley <- useDataset("hvulgare_eg_gene", mart = ensembl_plants)
zea <- useDataset("zmays_eg_gene", mart = ensembl_plants)

barley_gene_id <- "HORVU.MOREX.r3.4HG0381830"

getLDS(filters = "ensembl_gene_id",
values = barley_gene_id,
attributes = c("ensembl_gene_id","description", "go_id","name_1006"),
mart = barley ,
attributesL = c("ensembl_gene_id","description", "go_id","name_1006"),
martL = zea)

#> Gene.stable.ID Gene.description GO.term.accession
#> 1 HORVU.MOREX.r3.4HG0381830 NA GO:0016788
#> 2 HORVU.MOREX.r3.4HG0381830 NA GO:0016788
#> GO.term.name Gene.stable.ID.1 Gene.description.1
#> 1 hydrolase activity, acting on ester bonds Zm00001eb395980 Zm00001e038320
#> 2 hydrolase activity, acting on ester bonds Zm00001eb395980 Zm00001e038320
#> GO.term.accession.1 GO.term.name.1
#> 1 GO:0016788 hydrolase activity, acting on ester bonds
#> 2 GO:0016787 hydrolase activity
```
Question: I am using BWA to align files. I have four directories: seqtk_1, seqtk_2, seqtk_3, seqtk_4. Within each of those directories I have 10 subdirectories: subsample_1, subsample_2, subsample_3, etc. Within each of those subdirectories I have 20 paired-end reads (so from 10 genomes).

I want to put all the files from the all the directories through a pipe and into an output directory (BWA). The structure of this directory is the same as described above. So I have 800 input files and 400 output files.

I have written a script (below):

echo "[info] creating filenames";

for filename in ./Mock_Run/seqtk_*/subsample_*/*_1.fq.gz;
do file=`echo $filename|sed 's/_1.fq.gz//'`;
filenopath=`basename $file`;
for i in $(seq 10 $END);
do echo subsample_$i;
for u in $(seq 4 $END);
do eval outpath=BWA/seqtk_$u/subsample_$i;
echo "[info] starting BWA alignment...";
bwa mem -v 0 combine_reference.fa.gz ${filenopath}_1.fq.gz ${filenopath}_2.fq.gz > ${outpath}/${filenopath}_BWA.sam;
echo "[info] converting sam file to bam file";
samtools view -bS ${outpath}/${filenopath}_BWA.sam > ${outpath}/${filenopath}_BWA.bam;
echo "[info]filtering unmapped reads....";
samtools view -h -f 4 ${outpath}/${filenopath}_BWA.bam > ${outpath}/${filenopath}_unmapped.bam;
echo "[info] filtering mapped reads...";
samtools view -h -F 4 ${outpath}/${filenopath}_BWA.bam > ${outpath}/${filenopath}_mapped.bam;
echo "[info] sorting files";
samtools sort -o ${outpath}/${filenopath}_mapped_sorted.bam ${outpath}/${filenopath}_mapped.bam ;
samtools sort -o ${outpath}/${filenopath}_unmapped_sorted.bam ${outpath}/${filenopath}_unmapped.bam;
echo "[info] finished...no error to report";
done;
done;
done

It loops through all the files (like I wanted) and puts them into the right output subdirectory (like I wanted). It all seems to work, except it continues to loop. Once it has gone through all the files, it then starts again.

Any help would be appreciated.

Answer: You should parse the seqtk folder and subsample folder directly from `$filename` - see bellow for a suggestion. Your second and third for loops are already contained in the first:

First loop:

for filename in ./Mock_Run/seqtk_*/subsample_*/*_1.fq.gz;

Second loop:

for i in $(seq 10 $END);
do echo subsample_$i;

You are already globbing `subsample_*` and `seqtk_*`, so there is no need for `for i ... do subsample_$i` and `for u do seqtk_$u`, is there?

Try something like (untested):

SEQTK=$(echo $filename | sed "s|\.\/Mock_Run\/\(seqtk.*\)\/subsample_.*|\1|")
echo $SEQTK

Same for subsample.
Question: Hi guys,

The genotype imputation need human genetic map file, whose format is like(in chr1):

position COMBINED_rate(cM/Mb) Genetic_Map(cM)
55550 0 0
568322 0 0
568527 0 0
721290 2.685807669 0.410292036939447
723819 2.8222713027 0.417429561063975
723891 2.9813105581 0.417644215424158
728242 2.9806151254 0.430612871834774
729948 3.0780969498 0.435864105231133
740857 3.075133293 0.46941073432447

There are 3 columns, the first column is physical position united in bp, which is easy to understand. The second and third column is recombination rate(cM/Mb) and genetic position(cM, centiMorgan). I'm confused about the values in column 2 and 3.

Let's take the 4th row as an example. I know 721290 means the 721290th position on chr1. So what's the 2.68 and 0.41 represent for? I looked up some webpages and learned that "A centiMorgan is a unit of genetic distance that represents a 1% probability of recombination during meiosis." But I still cannot figure out the what's the meaning of 0.41/2.68 at position 721290.

Can anyone explain it to me? Thanks!

Tao

Answer: Hi, Tao.

Let's take the 4th row as an example.
1. 721290 means the 721290th position on chr1 and also means 0.721290 Mb.

2. To find 2nd column you need not numbers from 1st and 3d columns, but results of the subtractions: in the 4th row, it will be (0.41...-0)/(0.721290-0.568527) which gives us 2.685807669. So, in each row in column 1 and 3 you need to subtract numbers from the previous one to find the number from column 2.

Hope that it'll help you

Igor
Question: If I made a fake fasta with 10x the material of a book, randomly cut and spliced, and use this to align a Book, could tophat2 reconstruct the book? Or it only work with with ATCG letters?

Just a question that came out in my mind today.

Answer: You would need to re-encode the book as ACGT. For example, 1 ASCII character is 8 bits, corresponding to 4 nucleotides if you use the simplest possible encoding (rather than trying to pack into 7 or 6.5 bits, or whatever). Thus for an ASCII-formatted text file of the book, the encoded book would be 4x as long, but the mapping would work fine.

You MIGHT be able to map to the raw book using some protein aligners, as those allow more symbols.

As Istvan said, though, you'd need to use an assembler to reconstruct the book, not an aligner.
Question: Hello clever community!

I need your advice.
I am working on a *de novo* plant genome assembly of ~400 Mb.
I have Chromium 10x data, which was assembled with supernova. I also have Illumina paired end reads. Now I have additional data of PacBio reads, 120x roughly. The genome is diploid and I am thinking about using Falcon.

What do you think should be the best strategy:

1. Assembling PacBio reads and then using a tool to integrate the two assemblies? Is there anything like this? Which tool would you use?

2. Using a tool that can assemble the genome from both the chromium and the PacBio reads? Is there anything like it?

3. Assembling the PacBio reads and using chromium 10x and the illumina for polishing? If I assemble with Falcon, what tool should I use for polishing?

4? Anything else that I am missing to get the best out of what I can get?

Thank you very much in advance!
Alex

Answer: Falcon is not a bad choice, an alternative might be Canu (if you have the computational resources for it)

1) MEDUSA (as well as QuickMerge) is one of those integrating assembly/scaffolding tools

3) Pilon, Arrow, and there will be others I guess

4) Canu, but with the same remark as https://www.biostars.org/u/22207/ for MaSuRCa
Question: I am working on a project comparing RNAseq quantification results between Illumina short-reads and Nanopore long-reads and I have a couple questions about comparing the quantification results from these two technologies. More specifically I need some help with figuring out how to normalize the data for the comparisons within samples and between samples. So far I have come up with the following plan:

1. Using CPM to compare gene/transcript expression within each sample sequenced with nanopore. For example, comparing if gene.X transcripts are more abundant than gene.Y transcripts within sample_1 sequenced with nanopore. Using CPM instead of TPM for nanopore seems like a good option since our nanopore runs do not have transcript length bias. Does this sound like a good strategy?

2. Using TPM to compare gene/transcript expression within each sample sequenced with illumina. For example, comparing if gene.X transcripts are more abundant than gene.Y transcripts within sample_1 sequenced with illumina. Using TPM instead of CPM for illumina seems like a good option since illumina has transcript length bias (a single long transcript will have more counts that a single short transcript). Does this sound like a good strategy?

3. Here is where I am having trouble coming up with a good normalization strategy. Comparing gene/transcript expression between the same sample sequenced with illumina and nanopore. e.g., performing a spearman correlation between gene expression in sample_1 sequenced with illumina and sample_1 sequenced with nanopore. I am not sure what would work here since Illumina has transcript length bias and nanopore does not. Do you have any suggestions?

Any help here will be greatly appreciated.

Best, Bernardo

Answer: I did some work on this area before, it definitely has a lot of challenges.
The biggest difference, is short reads RNAseq measurement are calculated considering the transcript length (normalized against it).
The long reads tools does not normalized to transcript length. The reads are the actual transcripts, if you are using directRNAseq.
Due to the science behinds the chemistry.
I believe there is a paper that attempt doing this,

https://academic.oup.com/nar/article/50/4/e19/6439677

They use EM based method, similar to kallisto or salmon for illumina.
I will report back on my work later when I have more time.

I worked with the great @rpolicastro at the time.
He must have some great thoughts.
Question: Since phylogenetic analysis largely depends on a quality of primary sequence data I usually curate alignments manually and exclude phylogenetically uninformative or misleading sites.

But do you know any tool that could tell me how good my alignment is before sending it to phylogenetic program? Are there any programs or methods for assessing the phylogenetic signal that comes from a given alignment?

Answer: The 'quality' of an alignment is somewhat ambiguous given that an alignment is an inference of homology; we may or may not have a good idea of which sites are homologous in any given sequence set. Alignment algorithms compute alignment scores by assigning certain values to matches, mismatches, insertions/deletions, and gap extensions. These scores are then used to evaluate whether or not an alignment is better than another by simply comparing scores. However, the scoring scheme is arbitrary.

If you are concerned with the quality of your hand-curated alignment (and you may not need to be - expert 'by eye' alignments are often considered acceptable!), I would use your aligner of choice (MAFFT, perhaps?) and estimate the score of your alignment and compare it to the alignment produced or refined by the program.

One other concern: excluding misleading sites is one thing (the program Gblocks will remove regions thought to be resulting from spurious alignment), but removing non-informative sites can impact your analysis. For model-based phylogenetic approaches, invariant or slowly-evolving sites are included in the model as either a proportion of invariable sites or as part of the gamma distribution modeling among-site rate heterogeneity. If you are only selecting variable sites, it may be inappropriate to concatenate them and apply a model or a single evolutionary history. If this is the case, I would recommend [SNAPP][1] by David Bryant and others which estimates species trees from SNP data while treating individual gene trees as nuisance parameters.

Hope this helps.

[1]: http://www.beast2.org/wiki/index.php/SNAPP
Question: I want to compare two bam files to see whether they are same, how many reads they overlapped, and how many reads are unique. bam1 file was got several years ago and then I convert it to bam1_r1_fq and bam1_r2_fq and then remap to the new reference genome and got new sam file. Then I sorted, added header etc and got the bam2 file. When I use

samtools view -f 64 -F 2304 bam1 | cut -f 1,3,4 | LC_ALL=C sort -t ' ' -k 1,1 > bam1.txt
samtools view -f 64 -F 2304 bam2 | cut -f 1,3,4 | LC_ALL=C sort -t ' ' -k 1,1 > bam2.txt
join -t ' ' -1 1 -2 1 bam1.txt bam2.txt > bam.comparaison

I got error like this :

join: bam1.txt:4077: is not sorted: HWI-D00222:175:H7C06ADXX:2:1101:1250:20054	MT	10679
join: bam2.txt:4077 is not sorted: HWI-D00222:175:H7C06ADXX:2:1101:1250:20054	MT	10679

Then I want to try "cmpbams" that @Pierre Lindenbaum wrote, and it showed

mkdir -p lib/com/github/samtools/htsjdk/2.9.1/ && curl -Lk -o "lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar" "http://central.maven.org/maven2/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar"
/bin/bash: curl: command not found
maven.mk:102: recipe for target 'lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar' failed
make: *** [lib/com/github/samtools/htsjdk/2.9.1/htsjdk-2.9.1.jar] Error 127

How can I fix one of this problem or both? Thanks ahead.

Answer: Since many aligners work in non-deterministic mode the files will not be 100% identical (more so now due to different reference genomes, differences that may be there due to actual alignment logic, if you used different aligners/different versions of the same aligner etc).

It also appears that you are missing the program `curl` for your OS. Install it using appropriate means.