cannin/biostars_qa huggingface dataset:
Question: I am trying to get gene name from the SNP ID, I came across some posts and came up with this, but I am getting an error

Here is my code

library(biomaRt)
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl")
getBM(attributes = c("external_gene_name"),
filters = "snp_filter", values = "rs3043732", mart = ensembl)

This is the error

Error in getBM(attributes = c("external_gene_name"), filters = "snp_filter", :
Invalid filters(s): snp_filter
Please use the function 'listFilters' to get valid filter names

I am thinking that the "snp_filter" must have changed in the new version, but I am not able to find the new filter that serves the same purpose. Also is the "external_gene_name" is what the gene name that corresponds to the SNP ID ?

Answer: rs3043732 is no longer in dbSNP. Take a look HERE.

To look up rs IDs in biomaRt, you need to do this:

require(biomaRt)

ensembl <- useMart("ENSEMBL_MART_SNP", dataset = "hsapiens_snp")

getBM(attributes=c(
"refsnp_id", "chr_name", "chrom_start", "chrom_end",
"allele", "mapweight", "validated", "allele_1", "minor_allele",
"minor_allele_freq", "minor_allele_count", "clinical_significance",
"synonym_name", "ensembl_gene_stable_id"),
filters="snp_filter", values="rs6025",
mart=ensembl, uniqueRows=TRUE)

refsnp_id chr_name chrom_start chrom_end allele mapweight
1 rs6025 1 169549811 169549811 C/T 1
validated
1 1000Genomes,Cited,ESP,ExAC,Frequency,gnomAD,HapMap,Phenotype_or_Disease,TOPMed
allele_1 minor_allele minor_allele_freq minor_allele_count
1 C TRUE 0.00599042 30
clinical_significance synonym_name
1 benign,pathogenic,drug response,risk factor 17284
ensembl_gene_stable_id
1 ENSG00000198734

Look up all available attributes with `listAttributes(ensembl)`

Kevin
Question: Hello all,

this is a strange question perhaps, but it's brought on by a reviewer to one of my papers. We have used `mafft` to generate multiple sequence alignment for a simple phylogenetic analysis. However, the reviewer suggested we use "method that takes secondary structure of 16S" into account. I take it they meant something like `prank` aligner that uses codon information to generate a more accurate multiple sequence alignment. However, I don't know of any such tool for 16S.

Can you recommend anything relevant? Would be grateful for any suggestions.

Answer: Opinions differ about this, but MAFFT - or any "regular" sequence alignment program - is generally not able to align 16S sequences properly unless they are from fairly related species. It is not a matter of codon alignment - something that understands sequence co-variation is needed to align 16S sequences properly.

[**SSU-ALIGN**][1] can search for 16S sequences and automatically align them based on their phylogenetic membership.

[1]: http://eddylab.org/software/ssu-align/
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: sed 's/ *$//g' in.fasta > out.fasta

will remove only spaces at the end of lines. To remove tab or space use:

sed 's/\s*$//g' in.fasta > out.fasta
Question: Hi,

I am using the program Kallisto (verson 0.42.2.1) to perform pseudoalignment to a transcriptome and obtain estimated counts for each transcript. The program has an option to perform multiple bootstraps during this procedure, and I recently ran a test on one of my samples using 0, 1, 5, 10, 50, and 100 bootstraps to see at which point the estimated counts become relatively stable.

Surprisingly, the number of bootstraps did not change the estimated counts contained in the final output file abundance.tsv. However, this program also outputs an abundance file for each individual bootstrap run (compressed into abundance.h5), and the counts contained in these files were different. So, each time the program runs, the estimated counts were slightly different, but the end output file (abundance.tsv) did not change, no matter how many bootstraps were run.

If anyone could help with understanding why the estimated counts for each bootstrap would be different, but with no effect on the ultimate estimated counts, it would be appreciated!

Thanks,
Sarah

Answer: The output in the .tsv file is the maximum likelihood estimate for the expression of transcripts. That process is completely deterministic

The bootstrap outputs are stored in the .h5 files, they are generated by resampling the original data and finding the expression for this new dataset. This is why you will get slightly different est_counts in the bootstraps, and the amount of variation gives you an indication of how reliable the initial point estimate is.
Question: I have few motifs whose location on a given DNA sequence has been identified using MATCH (from TRANSFAC). MATCH provides the location of these matrices in a numerical way as motif 1 236 (+) (for an example - motif 1 is present at location 236). I want to get a graphical output by feeding the length of DNA sequence and the the position of matrices so that its easy for presentation and comparing two different sequences with almost similar motifs

Thanks in advance for any suggestions

Answer: You might want to look at berrylogo - it will be a more accurate means of visually comparing two motifs than a SeqLogo because a shift in nucleotide rank won't appear radically different

https://github.com/leipzig/berrylogo
Question: Hi I am trying to put together a plink ped file but I couldn't find how to depict missing maternal and paternal ids. You can give an argument if you have no parent ids in your file but I only miss some. For sex and phenotype you have missing data codes but what about the parental ids?
To put it differently how exactly do you depict the family relationships in the ped file format for plink?
Thanks

Answer: I got it. You just give 0 to parents maternal id and paternal id. This simply means that the parents' parents are missing and as long as you supply the true parental ids for the kids PLINK takes care of the rest.
Question: I am currently in the process of trying to understand what each clonality inference does under the cover and while doing that I realize that both the SciClone paper and the R-library feature [a multi-panel plot][1] to show Coverage vs VAF (Variant Allele Frequency) nicely aligned with a VAF histogram at the top. I also saw that the same plot type was also [featured][2] in a recent review article. However, neither paper gives a rationale regarding the information context for this plot and I would like to better understand what is unique about the VAF/Coverage graph and how it adds to our understanding of clones in addition to what we learn from the VAF distribution.

I am potentially over-interpreting things here, but just wanted to make sure that I don't miss any of the information resources.

[1]: http://journals.plos.org/ploscompbiol/article/figure/image?size=large&id=info:doi/10.1371/journal.pcbi.1003665.g001
[2]: http://www.nature.com/nrg/journal/v15/n8/fig_tab/nrg3767_F4.html
[3]: http://www.nature.com/nrg/journal/v15/n8/full/nrg3767.html

Answer: The coverage information isn't used in the clustering of VAFs at all (except to establish a minimum threshold). The plot has a couple of nice features though:

1) It can help identify sequencing artifacts (which often present with extremely high coverage and outlier VAFs)

2) it gives some "jitter" to the plots, so that the individual points can be seen in a 2d plot (otherwise they'd all be on a line)

3) Since there will be more sampling error with lower coverages, it helps ground your interpretation of the confidence of the clustering. (higher coverage experiment = more easily separable clusters and higher confidence in the results)
Question: I recently used PLINK to remove missingness in my SNPs data. I always saw message from terminal

Web-based version check ( --noweb to skip ) Recent cached web-check found...Problem connecting to web

Writing this text to log file [ HapMap_3_r3_4.log ] Analysis started: Thu May 6 20:03:31 2021

Although I used --noweb to skip connecting web, I am afraid I will use the web-based function in the future. Does anyone know if this is always happen in the PLINK? I'm using pLINK

Answer: This no longer happens if you upgrade to PLINK 1.9.
Question: Counting reads with htseq-count produces the following error

> Error occured when processing input (record #977823 in file
> <_io.TextIOWrapper name='' mode='r' encoding='utf-8'>):
> header not available in closed files [Exception type: ValueError,
> raised in libcalignmentfile.pyx:1884]

record #977823 is the last record (see below)

The command used to run htseq was

samtools view -F 2304 -f 1 -bh Samp1_PE.bam | samtools sort -n -O SAM - | htseq-count -s reverse -o Samp1_w_feature_PE.bam -p bam - /Annotation/Genes/masked.gtf >! htseq_pe.txt

The resultant text file is empty. The resultant bam file (Samp1_w_feature_PE) has all reads except for the last one (sorted by name). When counting the number of lines

> samtools view -F 2304 -f 1 -bh Samp1_PE.bam | wc -l

is exactly 977823, so I understand that record #977823 is the **last** record.

I have not found this error message when googling.
Any help?

Answer: Turns out that the specific parameter which causes this problem is

-o Samp1_w_feature_PE.bam -p bam

Counting paired-end reads AND creating the output in bam format causes an error.
None of them separately causes an error. I was able to output bam, when counting single reads, and more importantly, **I was able to output to sam format when counting paired end reads**. I have converted the sam to bam at a later stage, so this is effectively a successful workaround.
Question: Hi all, I need your help. I am struggling on getting a bar plot with ggplot2 package.

Imagine I have 3 different variables (which would be my y values in aes) that I want to plot for each of my samples (x aes):

str(df)
'data.frame':	4 obs. of 4 variables:
$ samples: Factor w/ 4 levels "S1","S2","S3",..: 1 2 3 4
$ v1 : int 234 189 167 235
$ v2 : int 345 265 235 254
$ v3 : int 54 48 32 645

I would like to have a barplot in which I have a group per sample (4), each of them with 3 bars, one for each variable (12 bars in total).
I know that if I want to extract the height from the df I have to call for stat = "identity".
So far this has been my guess, but it's not working properly:

n <- ggplot(df) +
geom_bar(mapping = aes(x = samples, y = v1), stat = "identity", position = "dodge") +
geom_bar(mapping = aes(x = samples, y = v2), stat = "identity", position = "dodge") +
geom_bar(mapping = aes(x = samples, y = v3), stat = "identity", position = "dodge")

Could you please help me?

Thanks before hand

Answer: with iris data : 4 observation of 4 variables:

data(iris)
library(ggplot2)
library(tidyr)

iris %>%
gather("Type", "Value",-Species) %>%
ggplot(aes(Species, Value, fill = Type)) +
geom_bar(position = "dodge", stat = "identity") +
theme_bw()

Separate boxes for each species:

iris %>%
gather("Type", "Value",-Species) %>%
ggplot(aes(Species, Value, fill = Type)) +
geom_col(position = "dodge") +
theme_bw()+
facet_wrap(~Species,scales = "free_x")

input (for one species, likewise for rest two species):

> head(iris)
Sepal.Length Sepal.Width Petal.Length Petal.Width Species
1 5.1 3.5 1.4 0.2 setosa
2 4.9 3.0 1.4 0.2 setosa
3 4.7 3.2 1.3 0.2 setosa
4 4.6 3.1 1.5 0.2 setosa
5 5.0 3.6 1.4 0.2 setosa
6 5.4 3.9 1.7 0.4 setosa
Question: Hi all,

I would like to download all protein sequences from one species on NCBI:

https://www.ncbi.nlm.nih.gov/protein?linkname=bioproject_protein&from_uid=261773

This is maybe trivial, but is there a way to download all sequences concatenated in only one fasta?

Thanks a lot,

Guillaume

Answer: using the ncbi interface you can just click on "Send to > File"

or using eutils:

```
curl "http://www.ncbi.nlm.nih.gov/entrez/eutils/elink.fcgi?dbfrom=bioproject&id=261773&linkname=bioproject_protein" | xmllint --xpath '//LinkSetDb' - | xmllint --format - | grep "" | cut -d '>' -f 2 | cut -d '<' -f 1 | while read L; do curl "http://www.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?db=protein&id=${L}&rettype=fasta" ; done
```

```
>gi|821074095|gb|KKY28990.1| putative uracil phosphoribosyltransferase [Diplodia seriata]
MFVHASGPESIKFKHLQGQVQVLLVDSVINSGATILDFVEAIREINPGIRIVVVAGTVQAQCISPNNPFY
KTLAQHGDISLVALRSSETKFTGSGGTDTGNRLFNTTHLL

>gi|821074094|gb|KKY28989.1| putative integral membrane protein [Diplodia seriata]
MPQYFPWPYSVDPLPEDLRRGLWPVGIFALMSTVATLALLCWITYRLVSWRKHYRSYVGYNQYVLLIYNL
LLADLQQSISFLISFHWIHTDSMLAPSPACFGQAWLVQIGDISSGMFVLAIALHTFFSVVKGRQIPFRAF
LIGTIVIWALALLLTVLGPALHGSDYFTAAGAWCWASDKYETERLWLHYLWIFIIEFGTVIIYALIFIYL
RKQLVSIASAHQHSTQNKVSQAARYMVLYPLTYVLLTLPLAAGRMATMTGQTLPIAYYCAAGSMMTSCGW
VDAALYALTRRVLVSNEIDQPQGGAGKGASSSGGRTGYGGHGSSHTATGWDIASFSDRKGGMGADHSVTI
TGGLDARGSNFIDMDELSKGGVHHHATERVGRPKHKGSSTPSTQGLTRARSSSTSARESTPRGSTDSILA
GLGGVRAETKVEIRVEPANGFMLPGEGSGSNGSSGMSTPNGRTVEVVGNSHAMRPRSGSPY

```
Question: Dear friends,

I am a cookie in RNA-Seq analysis so I am really confused with the Ensembl id and gene symbol. I have checked some associated posts in Biostars but I didn't find what I want eventually. The problem is that I noticed one Ensembl id matches multiple gene symbols, but I do not know how to deal with this issue when I want to do some analysis based on gene symbols. Should I add the counts of the same symbol together? Thanks in advance for the reply!

Answer: You have run into the problem that in the human genome there are instances of gene_names which are associated with multiple genomic loci (RF0019 in the link you posted). Since they are associated with different loci they also have different gene_ids. Last time I checked there were ~100 such gene_names in the human genome - many of which are located on different chromosomes.

I would always analyze the data with gene_ids (!) simply because else you assume the different loci produce identical products which might or might not be the case. Furthermore the gene_id analysis lets you analyze different things such as regulation and isoform switches. Lastly if you want to do any downstream analysis (go-terms or gene-set enrichment analysis etc) you should **NEVER** use gene_names. The problem is that in many cases gene_names are to unspecific with many different gene names pointing to the same gene and multiple genes all pointed to by a single gene name.

Cheers
Kristoffer
Question: I am trying to get gene symbols for gene ids that I got for mouse datasets. Gene ids look like that: `0610009B22Rik`. The code that I am trying to utilize is the following one:

ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
mouse_gene_ids <- dataset[, 1]
foo <- getBM(attributes=c('ensembl_gene_id',
'external_gene_name'),
filters = 'genedb',
values = mouse_gene_ids,
mart = ensembl)

I am getting `zero` results as an output after the query runs. I guess `filters` parameter is wrong. Any suggestions would be greatly appreciated.

Answer: The filter you need is `mgi_symbol` e.g.

library(biomaRt)

ensembl <- useMart("ensembl", dataset="mmusculus_gene_ensembl")
mouse_gene_ids <- "0610009B22Rik"

foo <- getBM(attributes=c('ensembl_gene_id',
'external_gene_name'),
filters = 'mgi_symbol',
values = mouse_gene_ids,
mart = ensembl)

Here's the result:

> foo
ensembl_gene_id external_gene_name
1 ENSMUSG00000007777 0610009B22Rik

I find the best way to choose the correct filter is to start with the Ensembl BioMart web interface, use the examples in the *Filters -> external references ID list* dropdown list to find the format I'm using, and then hit the XML button near the top. This will let you see the filter name required by **biomaRt**
Question: Hi All,

I went through [this paper on promoter atlas in FANTOM5][1], my question is that do they have cell wise CAGE dataset or is it global , as I cannot see the cell wise CAGE expression on their website.

Thanks in advance.

[1]: http://www.nature.com/nature/journal/v507/n7493/full/nature13182.html#the-fantom5-promoter-atlas

Answer: Hi:

Are you asking if they have cell-line/cell-type specific data to view? The answer is yes. A link to their ZENBU genome browser can be found here; you can add/remove certain datasets as tracks. While ZENBU works fine, for my use cases I prefer to use IGV for various reasons, including its flexibility. For example, you can load datasets on custom servers (ie iPlant Atmosphere) and it has integration with bedtools as well.

Hope this helps,

Taylor
Question: I have 10 folders, each containing a bam file and a corresponding bai file. I need to merge the bam files, which I believe I can do with `samtools merge merged_filename.bam *.bam`.

My question is that do I need to also merge the corresponding bai files as well? And if so, should I just do - `samtools merge merged_filename.bai *.bai`?

Answer: do

samtools index merged_filename.bam
Question: The recent UniProt website is fabulous, containing a complete collection of links to different resources. UniProt itself could be a very good resource for data mining for biological information.

I have a few questions regarding the SNP information in UniProt.

1. SNP information could be in two different section: Pathology and Biotech section and Sequence section. Is this from the annotation of disease SNP in dbSNP?
2. I want to identify disease SNP from natural variance. From xml file of UniProt (like: http://www.uniprot.org/uniprot/Q96CV9.xml), are all "feature" entries of variant with "evident" key disease related? Which is like following line: ``
3. How does UniProt do the mapping between dbSNP and UniProt? Just by sequence alignment of mRNA? Is the mapping complete and correct? Or there is any heuristic behind.

Thanks

Answer: 1. The disease kind is explained in http://onlinelibrary.wiley.com/doi/10.1002/humu.22594/full you can read a bit more on the other kind [here][1]
2. This is a bit easier using RDF and sparql http://www.uniprot.org/uniprot/Q96CV9.rdf than in the XML format. See [this example][2] to get all natural variants with a xref and involved in disease entry.
3. Hope one of my other colleagues comes round to answer that ;)

[1]: http://www.uniprot.org/help/variant
[2]: http://beta.sparql.uniprot.org/sparql?format=html&query=PREFIX+up%3A%3Chttp%3A%2F%2Fpurl.uniprot.org%2Fcore%2F%3E+%0D%0APREFIX+taxon%3A%3Chttp%3A%2F%2Fpurl.uniprot.org%2Ftaxonomy%2F%3E+%0D%0APREFIX+rdfs%3A%3Chttp%3A%2F%2Fwww.w3.org%2F2000%2F01%2Frdf-schema%23%3E+%0D%0APREFIX+skos%3A%3Chttp%3A%2F%2Fwww.w3.org%2F2004%2F02%2Fskos%2Fcore%23%3E+%0D%0ASELECT+*%0D%0AWHERE%0D%0A%7B%0D%0A%09%3Fprotein+a+up%3AProtein+%3B%0D%0A%09%09up%3Aorganism+taxon%3A9606+%3B%0D%0A%09%09up%3AencodedBy+%3Fgene+%3B%0D%0A%09%09up%3Aannotation+%3FdiseaseAnnotation+%2C%0D%0A++%09%09%09%3FnatvarAnnotation+.%0D%0A++%09%3Fgene+skos%3AprefLabel+%3Fname+.%0D%0A%09%3FdiseaseAnnotation+a+up%3ADisease_Annotation+%3B%0D%0A%09%09rdfs%3Acomment+%3Ftext+.%0D%0A++++%3FnatvarAnnotation+a+up%3ANatural_Variant_Annotation+%3B%0D%0A++++++++rdfs%3AseeAlso+%3Fdbxref+.%0D%0A%7D
Question: Hi,

I am trying to convert a vcf.gz file to bed/bim/fam format using PLINK2 (data is hg38).

My command is:

plink2 --vcf variants.vcf.gz --out variants --freq --make-bed --chr 1-22 XY

I get an error:

Error: Invalid chromosome code 'chr11_KI270721v1_random' on line 36088290 of

How can I ensure that I filter out these variants/sequences?

thanks!

Answer: Add the "--allow-extra-chr" (or "--aec" for short) flag to permit these codes in the input. (--chr 1-22 XY will then filter them out.)
Question: Hello, everybody!

I am looking for the procedure reversed to cat.

I would like to separate a long fasta-file into many single fasta sequences, so each sequence from the large fa-file will be situated in a new small file.fa each. I usually used something like the following, but it not as good as it looks like... Where is the mistake?

Many thanks for your help!

N.

```
sub get_fasta {
my $inputName = $_[0];
my $fh;
# new input record separator for fasta-files
$/ = "\n>";
open ($fh, $inputName) or die "cannot open file: $!\n";
my @entries = <$fh>;
# remove last character '>' except for the last entry
chop(@entries[0..$#entries - 1]);
# add '>' except for the 1st entry
foreach (@entries[1..$#entries]) {$_ = '>'.$_;}
# print"@entries";
close($fh);
$/="\n";
return \@entries;
```

Answer: Another option in bash... Possibly more concise and readable...

while read line
do
if [[ ${line:0:1} == '>' ]]
then
outfile=${line#>}.fa
echo $line > $outfile
else
echo $line >> $outfile
fi
done < myseq.fa

You input file is myseq.fa. Individual files will have the name of the corresponding sequence, without leading `>`. If your sequence names contain funny characters (non alphanumeric) it will probably break but could be fixed depending exactly how you want to call the output files.
Question: I would like to get a nonredundant set of heterodimers (30% seq.identity) with PDB ids and its corresponding chains. Is it possible to get with pdbePISA? I tried , but I couldn't find. Your suggestions would be appreciated.

Answer: Try the new stoichiometry drilldown at www.rcsb.org . Click on 99293 at the top of the page, then "Protein Stoichiometry"-> Heteromer->Hetero-2-mer , then "Show only representatives at 30%"
Question: I have two files:

1. An annotation file in gff format
2. A fasta formatted sequnce file for a whole genome.

I want to have a resulting file : myGenomeCDS.fasta

How can I go about doing this via python or Tools ?

Thanks

Answer: There are so many ways to do this. If you don't want to write your own code, use the utilities from preexisiting packages such as BedTools (fastafrombed), Glimmer (extract) and even Galaxy as a tool (extract features).
Question: I am trying to calculate the forward and reverse mapped reads from bam. I am eventually going to use a bed of targeted regions but started by using a script from a thread from @Alex Reynolds.

samtools view Reads.bam | gawk '(and(16, $2))' > forwardStrandReads.sam
samtools view Reads.bam | gawk '(! and(16, $2))' > reverseStrandReads.sam

output for forwardStrandReads

1QXCL:01920:06690	16	chr1	10177	0	34S41M	*	0	0	TCGTTACCTACCTACCTACCTACCTAACCTACCTACCTAACCCTAACCCTAACCCTAACCCTAACCCTAACCCTA	+6+(+.*44456775:44-:44-74/5/:773855/9965*78737*88848+9994<4?@>8<6;<<7808987	ZP:B:f,0.00718598,0.00485627,0.000543075	ZM:B:s,286,-24,286,12,14,256,-6,274,10,206,226,0,488,514,18,2,250,0,24,494,16,268,222,238,498,0,2,234,4,-26,658,-18,468,218,-28,700,476,260,6,734,456,6,-6,-14,242,682,28,-30,454,-2,-24,250,8,16,658,18,440,2,-12,236,34,-4,648,0,472,252,-12,564,346,210,48,574,332,10,-24,-12,244,530,56,-10,398,-34,30,224,-26,84,544,20,344,42,-38,228,82,-22,524,-26,340,222,-30,484,282,192,56,526,236,46,-22,0,320,482,76,12,324,14,48,382,80,180,258,168,124,102,-70,96,310,54,314,10,178,226,54,224,110,126,160	ZF:i:24	MD:Z:41	NM:i:0	AS:i:41	XM:i:41	XA:Z:map4-1	XS:i:41	RG:Z:1QXCL.IonXpress_001-7EBD8685	PG:Z:tmap-1F1907BA

I guess my question is what am I looking for as I am nor sure how to interpret this. Is there documentation that might help? Thank you :).

Answer: Hi cmccabe,

just go to [http://broadinstitute.github.io/picard/explain-flags.html][1]. There you enter the number seen in the second column. In the given case, that's 16. It'll give you a tick at "read reverse strand".

If you are dealing with RNA-Seq data, have a look at [RSeQC's infer_experiment tool][2], since it copes with paired-end information and the read-orientation with regard to the transcripts' orientation.
Otherwise, have a look at [bam_stat][3], it gives you the number of reads mapping to each strand.

Cheers,

Michael

[1]: http://broadinstitute.github.io/picard/explain-flags.html
[2]: http://rseqc.sourceforge.net/#infer-experiment-py
[3]: http://rseqc.sourceforge.net/#bam-stat-py