cannin/biostars_qa huggingface dataset:
Question: I am new to the protein structure prediction, so the following question may sound obvious to experts in the field.

Let's say a protein has 3d crystal structure already, and alphafold is used to predict the protein structure and generates a .pdb file. How do we know how close the predicted 3d structure is to the 3d crystal structure?

Answer: The relationship of the two structures is tested by [**structural alignments**][1]. There are many servers and programs to do that.

I should tell you that this is not really a way to test the accuracy of AlphaFold, or any prediction program for that matter. Unless you excluded the identical template - and all related templates - chances are that AlphaFold used it for modelling, and the prediction should be very accurate. In terms of structural alignments, that means low [**root mean square deviation**][2].

Others who work in that field have already tested AlphaFold in an unbiased fashion (yours would be decidedly biased), and found it to be very accurate. The program itself will give confidence values for its models. Rather than doing what you described, I suggest you try it on a protein without a known 3D structure.

[1]: https://en.wikipedia.org/wiki/Structural_alignment
[2]: https://en.wikipedia.org/wiki/Root-mean-square_deviation
Question: In GWAS studies, 5*10-8 seems the defacto standard for association. A lesser known threshhold is 1*10-5.

This threshhold is often called "suggestive association."

I am writing to ask, where did this figure come from? Were empirical tests or statistical simulations done to select this level?

I can find many papers that USE this level, but far fewer that JUSTIFY it. Any relevant citations would be welcome, and thank you.

Answer: > "The P values of Bonferroni corrected thresholds for suggestive, 5 and
> 1% genome-wide significant levels were 1, 0.05 and 0.01, respectively,
> divided by the number of SNPs used in the GWAS. The suggestive level
> was first proposed by Lander and Kruglyak [17] and represents the
> threshold where, under the null hypothesis, one false positive is
> expected per genome scan."

So the suggestive line is calculated as the `-log10( 1 / number of variants) `

[Reference](https://www.ncbi.nlm.nih.gov/pubmed/28196480): Guo, Y., Huang, Y., Hou, L., Ma, J., Chen, C., Ai, H., ... & Ren, J. (2017). Genome-wide detection of genetic markers associated with growth and fatness in four pig populations using four approaches. Genetics Selection Evolution, 49(1), 21.
Question: Hello Biostars,

I have a question that is perhaps naive. I use hclust() to cluster my data and heatmap.2() to map, then use cutree() to get subclusters. But the order of subclusters I got from cutree() is not the same as the order visualized on the map. I know the problem is related to the order. I have read some materials but still feel confused. How to solve it?

Thank you!

Answer: To get the subclusters of your data after clustering, use the vector returned by cutree as an index into your original data matrix. Here is a simple example:

```
library(gplots)

# create some data
d <- matrix(rnorm(120),12,10)

# cluster it
hr <- hclust(as.dist(1-cor(t(d), method="pearson")), method="complete")

# define some clusters
mycl <- cutree(hr, h=max(hr$height/1.5))

# get a color palette equal to the number of clusters
clusterCols <- rainbow(length(unique(mycl)))

# create vector of colors for side bar
myClusterSideBar <- clusterCols[mycl]

# choose a color palette for the heat map
myheatcol <- rev(redgreen(75))

# draw the heat map
heatmap.2(d, main="Hierarchical Cluster", Rowv=as.dendrogram(hr), Colv=NA, dendrogram="row", scale="row", col=myheatcol, density.info="none", trace="none", RowSideColors= myClusterSideBar)

# cutree returns a vector of cluster membership
# in the order of the original data rows
# examine it
mycl

# examine the cluster membership by it's order
# in the heatmap
mycl[hr$order]

# grab a cluster
cluster1 <- d[mycl == 1,]

# or simply add the cluster ID to your data
foo <- cbind(d, clusterID=mycl)

# examine the data with cluster ids attached, and ordered like the heat map
foo[hr$order,]
```

If you run the code, you'll see that mycl returns the cluster membership for each row of your original data. To see it in the order as it is plotted in the heatmap, you have order it according to the heat map index. Another caveat many find confusing, is that the first position of the vector refers to the bottom of the heat map.
Question: I have 10 gvcfs for 10 samples and a list of SNP positions (chr# and position) I want to extract from the gvcfs. I can't just use vcftools to do this because the position is given as a range (start position and end position) in the gvcf and if the position in my list is within the range but not the start of the range, it won't be extracted. How can I combine my 10 gvcfs into a single vcf with only the SNPs at the positions I specify in my list?

Answer: you don't do that. from a gvcf you call GenotypeGVCFs https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_tools_walkers_variantutils_GenotypeGVCFs.php using the option `-L` for the desired position(s). ( https://software.broadinstitute.org/gatk/documentation/tooldocs/3.8-0/org_broadinstitute_gatk_engine_CommandLineGATK.php )
Question: Hello all,

I'm beginning to work with data bases form my Master Thesis.

I want to map each organism name with the corresponding UCSC Genome Browser assembly ID to save it in a table. As I don't have the UCSC database available to work with MySQL my question is if there is any way to get all the UCSC Genome Browser assembly ID.

Thanks you all!

Answer: Genomax already answered, but here is also a mysql command you can use if you don't want to download the entire hgcentral database dump:

$ mysql -h genome-mysql.soe.ucsc.edu -u genome -e "select name, scientificName from dbDb limit 5" hgcentral
+---------+-------------------------------------+
| name | scientificName |
+---------+-------------------------------------+
| hg38 | Homo sapiens |
| balAcu1 | Balaenoptera acutorostrata scammoni |
| calMil1 | Callorhinchus milii |
| rn6 | Rattus norvegicus |
| sorAra2 | Sorex araneus |
+---------+-------------------------------------+

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

- General questions: genome@soe.ucsc.edu
- Questions involving private data: genome-www@soe.ucsc.edu
- Questions involving mirror sites: genome-mirror@ose.ucsc.edu

ChrisL from the UCSC Genome Browser
Question: Hi all,

I've a simple but intriguing question: I was trying to pull the BLOSUM80 matrix online, and I see that the substitution scores differ among sources.

For example, the Ter scores on [this matrix][1] are `-8` while the same change on this [NCBI BLOSUM80][2] are `-6`

Is there a "right" matrix? What do these different scores indicate as an underlying factor? I'd appreciate the community's input on this please.

Thank you!

[1]: http://www.cbcb.umd.edu/confcour/Fall2008/CMSC423-materials/BLOSUM80
[2]: http://www.ncbi.nlm.nih.gov/IEB/ToolBox/C_DOC/lxr/source/data/BLOSUM80

Answer: It looks like the two versions are using different units. In the first matrix the unit is 1/3 bit, while in the second matrix, it's 1/2 bit.

From the second matrix:

a scale of ln(2)/2.0

is the conversion factor for 'nats' to half bit units.
Question: I have some short reads aligned using BLAT, the output is in tabular psl format (including the sequences) for each alignment. Is it possible to convert the blat output to SAM/BAM format. Myself, I would think it is not because of the lack of some data fields in psl which is required for SAM format (mainly the CIGAR string), but please proove me wrong!
Normally I would advise myself to use a different tool (bwa, bowtie, lastz) and align getting a SAM file, but what if that wasn't an option (say because you really want to use BLAT or you don't have the input) is there a way to do the conversion. I can possibly code that in perl and share it if someone had an idea how to do it.

Answer: (not tested): there is a script named `psl2sam.pl` in:

samtools/misc/psl2sam.pl
Question: ~~I am trying to update my bioconductor to the latest version (3.0+). I currently have R 3.2.1 and Bioconductor 2.28. When I run BiocUpdate, it just reinstalls 2.28. I tried uninstalling and reinstalling Bioconductor, but still only got 2.28.~~

~~The reasoning is that a package I am running (phyloseq) is having an error that the developer doesn't recognize. He suggested updating to the latest version of Bioconductor & R to fix it. I have updated my R version, but am unable to update my Bioconductor.~~

~~Any suggestions to install Bioconductor 3.0 or 3.1?~~

**Edit: There was a misunderstanding in how to tell what version of Bioconductor I have.**

**The best way to check the version of Bioconductor is to compare your version of BiocInstaller (a package within R once Bioconductor is installed) with the Bioconductor website.**

Answer: As of about R 3.2.0, the way to get the associated Bioconductor version is to run the non-exported built-in function:

```r
> tools:::.BioC_version_associated_with_R_version()
[1] '3.1'
```

Mind the `.` in front of the `BioC`.

That said, it may be possible for some installed packages to have been sourced from a different BioC version.
Question: I'm struggling to represent this syntax with CWL:

vcf-sort -c in_file.vcf > out_file.vcf

Specifically with routing the outfile, as right now, it looks like this will push to stdout.txt.
I'm sure this is simple, but, guidance is appreciated.

cwlVersion: v1.0
class: CommandLineTool
requirements:
- class: DockerRequirement
dockerPull: preprocess
inputs:
vcf:
type: File
inputBinding:
position: 1
prefix: -c
outputs:
- id: output
type: File
outputBinding:
glob: "sorted-output.vcf"
baseCommand: vcf-sort

Answer: outputs:
- id: output
type: stdout

stdout: sorted-output.vcf

:-)
Question: What is the easiest method in 2017 to obtain build 37 positions for a batch of rsIDs? This site does not seem to support specific builds:
https://www.ncbi.nlm.nih.gov/projects/SNP/dbSNP.cgi?list=rslist

The opposite question has been asked in here; go from coordinate to rsID.

I wouldn't mind Python based suggestions at all.

EDIT: I found a solution using the Entrez class of the Python module Biopython here:
http://www.danielecook.com/getting_snp_dat/

It's a bit convoluted. I'll simplify it a bit and post the solution.

Answer: ```
$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -D hg19 -e 'select chrom,chromStart,chromEnd,name from snp147 where name in ("rs371194064","rs779258992","rs26","rs25")'
+-------+------------+----------+-------------+
| chrom | chromStart | chromEnd | name |
+-------+------------+----------+-------------+
| chr7 | 11584141 | 11584142 | rs25 |
| chr7 | 11583470 | 11583471 | rs26 |
| chr1 | 10149 | 10150 | rs371194064 |
| chr1 | 10146 | 10147 | rs779258992 |
+-------+------------+----------+-------------+
```
Question: Hello biostars,

I have pre-processed mRNA data from an illumina whole-genome expression array Human HT-12 V4 and microRNA data from an Agilent GeneSpring GX 11 array.

I am aware that, low level expression genes (<10) should be removed from RNA-seq analysis.
Is there such a threshold or requirement to do so with microarray data.

I have not come across any literature that says so, but I would like to know if others remove low level genes from microarray data.

Any help or advice is appreciated.

Answer: In general, microarray probes should be filtered if their intensity is not above the background intensity - how to detect this varies by array platform.

Check the limma [User Guide][1], it shows how to filter "non-expressed" probes in several cases. For example, for the Illumina BeadChip Arrays - like the Human HT-12 V4 - the User Guide says:

>Filter out probes that are not expressed. We keep probes that are expressed in at least three arrays according to a detectionp-values of 5%:
>
> expressed <- rowSums(y$other$Detection < 0.05) >= 3
> y <- y[expressed,]

See page 111 of the User Guide (revision of 23 January 2019), read both **17.3.4 How many probes are truly expressed?** and **17.3.5 Normalization and filtering**.

If you read this guide carefully, it will show how to filter noisy probes for other types of arrays as well.

[1]: http://www.bioconductor.org/packages/devel/bioc/vignettes/limma/inst/doc/usersguide.pdf
Question: I have a bed file with transcript_chr/transcript_start/transcript_end/....../snp_pos/snp_type fields. Multiple SNPs may belong to the same transcript (so multiple rows may have the same trarnscript_start/transcript_end). I would like to use python or R to plot the transcripts maybe as rectangles and SNPs as points over the rectangles to see the distribution. Is there a specific tool (package) i could use for that (R or python)?

Thank You!

Answer: GenomeGraphs might be helpful in R, check Fig 2 here

https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-10-2
Question: I am working on a project to detect positive selection in several different genes.

I am working with colobus, rhinopithecus, macaca, gorilla, pan, pongo, homo, andrillus, papio, and have dogs and mice in there as well. We plan to separate these organisms into foreground and background clades as per Yang et al 2005 (Evaulation of an improved branch-site likelihood method for detecting positive selection at the molecular level). I am then going to attempt to detect positive selection using the branch site model test described in that same manuscript.

I feel comfortable with the statistical and theoretical parts of the exercise, however, I could really benefit from expert guidance on selecting sequences. Specifically, I am wanting to ask:

1. Beyond resources like genbank, is there a good place to get sequences on several different genes in each of these species? When doing this, is it acceptable retrieve sequences from a resource like that if you intend on publishing a positive selection paper, or is there a recommended protocol? Or are there multiple databanks?
2. My second question revolves around making the correct phylogenetic tree for the above organisms. I have obtained this tree from Clustal, have sanity checked it, and it makes sense with what I understand to be true, but I think that in the literature there is a standard way of making these trees. Could anyone briefly outline the best way to construct a tree for these organisms? With any luck it will validate what I already made, but I definitely want to check that.
3. Last, I would like to call on the experience of those who have conducted tests for positive selection based on dN/dS statistics: what are some of the most common or sinister ways of generating false positive results? Anything to double check to make sure my sequences and alignments are solid so that the statistics are measuring positive selection, and not some artifact of the sequence selection or alignment protocol?

Thank you very much for your time and help.

Answer: 1. You can look at pulling down genes for sets of taxa using the [PhyLoTA browser][1]. Alternatively, you can try to get the genomes and extract the genes directly (not sure if you have genomes available for all of those). But, the best way is probably write some code to query Genbank and look for what you want, especially if you're focusing on generating a lot of single-locus datasets.
2. You want to thoroughly estimate a phylogenetic tree, not just use one that's probably being used as a guide tree for alignment in Clustal. This has been answered a lot on Biostars and other resources, but you can start with one of my older posts here: https://www.biostars.org/p/113346/
3. The power to detect positive selection using tree-based approaches scales with the total amount of divergence among taxa in the tree. Roughly, you have more power to detect selection if your tree is 'deep' as opposed to 'shallow.' Check out some of [Maria Anisimova's work][2] and most model/PAML papers by Ziheng Yang and others.

Let me know if you have any other questions.

[1]: http://phylota.net
[2]: http://mbe.oxfordjournals.org/content/19/6/950.long
Question: I have a bigWig file with the following profile:

bigWigInfo uniqueness.bw
version: 4
isCompressed: yes
isSwapped: 0
primaryDataSize: 919,783,047
primaryIndexSize: 97,088,256
zoomLevels: 10
chromCount: 25
basesCovered: 3,095,693,183
mean: 0.828284
min: 0.000000
max: 1.000000
std: 0.361063

I would like to create a bed for all regions with value < 1.

More generally, I'd like to be able to threshold a bigWig on a value and get all regions above or below that value in bed format.

Thanks a lot!

Answer: Copy-pasting from [where you cross-posted](https://bioinformatics.stackexchange.com/questions/2800/bigwig-to-bed-for-regions-above-below-threshold/2801#2801):

You can do that with a bit of python:

#!/usr/bin/env python
import pyBigWig
threshold = 10 # Change me
bw = pyBigWig.open("file.bw") # Change me
of = open("regions.bed", "w") # Change me

for chrom, len in bw.chroms().items():
intervals = bw.intervals(chrom)
for interval in interals:
if abs(interval[2]) > threshold:
of.write("{}\t{}\t{}\n".format(chrom, interval[0], interval[1]))
bw.close()
of.close()

You'll need to install pyBigWig (it's available via pip or conda install) and change the lines with Change me in them.
Question: Hi

Extracting only first sequence from the multiline fasta file

perl -02 -ne '/(>[^>]+)/ && print $1' OG0012884.fa

>TRINITY_DN9_c0_g1_i1.p1
MPMKGRFPIRRTLEFLRSGTVVFKDSVKILTVNYNTHGERSDGARKFVFFNIPQIQYQNP
WIQILMFKNMTPSPFLRFYLDDGEQVLVDVEGKNHKQIVEHVKTILGKNDVLLEADKQVQ
KEHSHPAHFGPKTYCLRECMCEVGGQVPCPGVVPLPKEMTGKYWTALRAGSAI*

How do I extract only the 2nd and 9th sequence, ?

Using Bash able to print only 1st ID in a fasta file ?

head -1 OG0012884.fa
>TRINITY_DN9_c0_g1_i1.p1

How do i extract only the 2nd ID ?

Help is need, suggestions please.

Answer: Here's an awk based solution which counts the headers as they come and prints the 2nd and 9th along with their sequence on the respective next lines.

awk '{if ($1 ~ "^>") count+=1} { if ( count == 2 || count == 9 ) { print $0; t = NR+1};} { if ( NR == t ) { print $0 }}' OG0012884.fa

If you happen to be working with a file that contains the sequence on the same line:

awk '{if ($1 ~ "^>") count+=1} { if ( count == 2 || count == 9 ) print $0}' OG0012884.fa

Heres the pseudocode, if you care to follow along.

for line in file:
if line starts with ">":
count += 1
if count == 2 or count == 3:
print line
t = current line number + 1
if current line number = t:
print line
Question: Dear all,

I have a molecule in PyMol which I want to open in ChimeraX. But after exporting the molecule as .pdb, some structural information seems to be lost - regions that in the PyMol session file are embedded into helices are now loops (this regards to several residues only). If I open the .pdb in PyMol, it looks the same as in chimeraX. I have no idea what I keep doing wrong, any hints? The first 14 residues of the original protein (from PDB, corresponding to Histag) were removed with a custom PyMol script which changed residues enumeration so that it properly starts with 1 - that is why the molecule from PyMol is needed. I would be grateful for the help.

best,

maciek

Answer: PDB files have information about secondary structure elements that looks like this:

HELIX 11 11 GLY A 219 ASP A 222 1 4
HELIX 12 12 GLN A 235 CYS A 240 1 6
HELIX 13 13 TYR A 248 SER A 252 1 5
SHEET 1 A 6 ASN A 53 HIS A 57 0
SHEET 2 A 6 VAL A 65 THR A 69 -1 N THR A 69 O ASN A 53
SHEET 3 A 6 VAL A 29 GLN A 33 -1 N LEU A 32 O ALA A 66

Most visualization programs will use this information when rendering molecules. However, I don't know if this information is preserved when you save the file in PyMol. If not, PyMol will use its internal algorithm to calculate helices and strands, which can give different results from the original file. I think something similar could be true for Chimera.

Long story short: I suggest you open the original file in a text editor and delete the HIS-tag by hand - you will have to scroll down a while to find `ATOM` residues and delete the histidines. As long as you don't change anything else, you should be able to get the original secondary structure elements when you open the file afterwards. If you don't know how to do this, open the original PDB file and find all the lines that begin with either `HELIX` or `SHEET`. Copy and paste them into the file created by PyMol, and make sure they are right before the line that starts with `CRYST1`.
Question: I'm trying to do a bedtools compare on a couple of bed files. One of the bed files keeps throwing this error:

Error: unable to open file or unable to determine types for file

My file is as follows:

```
chr1 810865 3198369
chr1 844270 845356
chr1 882432 10373009
chr1 1104962 1173985
chr1 2616309 2617058
chr1 3056425 3245459
chr1 4704545 5447621
```

I have checked that the start position is before the end position. I have added and removed column headings chrom, "start" and "end" (also tried 'chromStart' and 'chromEnd'. But no joy. I had originally created this bed file from a tab delimited file from R. Does anyone know what has gone wrong?

Answer: OK so this is what I did in case anyone else gets the same problem

1. As per Ashutosh I ran `perl -p -i -e 's/ /\t/g'` I think this was not the problem though in the end but is probably a good idea
2. Uploaded the bed file to UCSC genome browser as a custom track. If the bed file is incorrectly formatted it will tell you what the problem is and is much better at debugging than bedtools
3. beware of directory paths when running the command line.

Sorted now

Thanks everyone
Question: Hi

Can anyone help in SNP analysis. Which online tool should I use to get the data for SNP.

Answer: Owhkee, so you did something like create a bowtie database (or bwa or whatever mapper you used), then aligned your reads or query against your indexed reference. Then you piped that to a Sam or Bam file right?

Well now you can do SNP calling, there are a lot of tools, i use Freebayes but I have also used Samtools Mpileup and BCFtools to get SNPs but for me freebayes works the best.

If you have more specific questions just aks, but i would first just read a little about SNP calling tools so you know for sure that you have the right tool for your project. And read the manuals to understand what the tools do. For freebayes you can look at [this tutorial](http://clavius.bc.edu/~erik/CSHL-advanced-sequencing/freebayes-tutorial.html).