cannin/biostars_qa huggingface dataset:
Question: Hello everyone,
I would like to ask you how can I get the number of exon? I want visualize gene length vs. number of exons in gene.

For example:

I have the length of the gene:

exons <- getBM(attributes = c("external_gene_name","start_position","end_position","ensembl_exon_id"),
filters = "chromosome_name",
values = "21",
mart = genes)
exons[5] <- exons$end_position - exons$start_position
names(exons)[5] <- "gene_length"

I tried to find out the number of exons:

as.data.frame(table(exons$ensembl_exon_id))

But now, I don't know what else. Please, somebody would help me?

Answer: Something like :

require(tidyverse)
genes <- useDataset("hsapiens_gene_ensembl",mart=ensembl)
exons <- getBM(attributes =
c("external_gene_name","start_position","end_position","ensembl_exon_id"),
filters = "chromosome_name",
values = "21",
mart = genes)
exons <- as.tibble(exons)
exons.plot <- exons %>%
mutate(len=end_position-start_position) %>%
group_by(external_gene_name) %>%
summarize(geneLen=sum(len),exonNumber=n())
# plot
ggplot(exons.plot,aes(x=geneLen,y=exonNumber))+
geom_point(alpha=0.5)+
scale_y_log10()+
scale_x_log10()+
theme_bw()

FYI you should use ensembl_gene_id instead of gene name

![enter image description here][1]

[1]: https://image.ibb.co/byyCpb/Rplot02.png
Question: I am looking for a tool that can provide a high level summary of coverage statistics of genomes used in an alignment.

I concatenated GCF_000002285.5_Dog10K_Boxer_Tasha_genomic.fna, GCF_000001405.39_GRCh38.p13_genomic.fna, and GCF_000001635.27_GRCm39_genomic.fna. Each of these animal genomes contains hundreds of chromosome sequences.

I have aligned RNASeq data to this combined dog-human-mouse reference and now I would like to calculate the total number of reads, coverage of each animal genome and depth of coverage. I am only interested in animal level summary, not chromosome coverage statstics.

Tools such as mosdepth, depth-cover, qualimap (V2.2.1) only provide coverage at a chromosome level, not a summary at full genome level (eg coverage 1x, 4x, 10x depths).

>Update after a couple of days: Pehaps what I can do is get all the chromosome accessions for each of the genomes. Then extract from the .bam (forst convert to sam) result 3 seperate files containing only the results relevant to each genome. Then find a tool that will calculate summary statistics on each of the 3 sam files (total coverage at depth n, total number of reads mapped).

[1]: https://github.com/shiquan/bamdst
[2]: http://bamstats.sourceforge.net/

Answer: Will use [Pysam][1] to write code to do this

[1]: https://pysam.readthedocs.io/en/latest/api.html#api
Question: Hi there!

I'm looking only for the records of **nuclear** loci for a certain plant species within `Nucleotide database` on the `NCBI` website.Options to select either chloroplast or mitochondrion or plastid are available. I'm wondering, if it's would be possible to easily select only **nuclear** loci on the website as well?

Many thanks,

Denis

Answer: You can build a query to exclude the genetic compartments such as chloroplast, mitochondrion, and plastid. `"viridiplantae"[Organism] NOT "plasmid"[Filter] NOT "plastid"[Filter] NOT "mitochondrion"[Filter]`
Question: Hi all,

I'm having problems with RNA-seq data. I have two replicates for five different treatments (all on the same cell line) and also two replicates for the control. I used DESeq2 to perform the differential expression analysis and I obtain that all the treatments are very similar to each other (I mean that the percentage of differentially expressed genes is always around 1 to 3% for each treatment). Only for one treatment I found that the percentage of differentially expressed genes is around 20%. This was unexpected and I can't find out a method to evaluate if it is a biological or a statistical issue. Do you have any idea on which analysis I could perform to confirm that I carried out correctly the differential expression analysis? If it is a biological issue do you have any idea of a possible explanation of such a substantial difference?

Thanks,
AriBo

Answer: > ...percentage of differentially expressed genes is always around 1 to
> 3% for each treatment). Only for one treatment I found that the
> percentage of differentially expressed genes is around 20%. This was
> unexpected...

What was unexpected about it? These numbers are not inconsistent with what we know is possible biologically. What is it that you know about your experiment, that leads you to believe this is an odd result? Given the information in your question there's nothing incompatible with a biological explanation. The PCA plot sets these samples apart from the others, so a large DE response could make sense. If you look at the values for DE genes and non-DE between each condition and control, does the data makes sense to you? Is there a core set of controls or housekeeping genes that exhibit similar behavior/characteristics across the experiment?

Is there anything about the biology in your results that does make sense? Or does NOT make sense (i.e. you applied heat shock, but your cells are exhibiting a cold shock response)? Drawing a line between each condition and at least one or a few genes one might expect to respond should be doable. Is there a confounding factor you haven't mentioned? For instance, was everything controlled for except condition? (i.e. all cells and conditions prepared and applied on the same day in the same laboratory by the same person?) or is there a batch effect whereby the outlier condition was applied 3 months later in a different laboratory on the other side of the world (but maybe nobody mentioned that)?
Question: Hello,

I have created a local blast database and blasted some genes against the reference using

blastn -query a.fasta -db BLASTDATABASE -max_target_seqs 5 -outfmt 6 -out output

I obtain a ~98% similarity score for the best hit. Now, I want to extract the `fasta` sequence of that BLAST search together with the 200 bp flanking sequence from both start and stop sides for designing a primer.

Can you please give me a suggestion on how to do this for a single BLAST search?

Thank you.

Answer: For extracting the region + context of the hit, you're best of with using `blastdbcmd`. it has an option `-range` that you can use to extract a certain region from a hit in the blastdb .

As with the approach https://www.biostars.org/u/50673/ proposed (for the query) you can extract the start & stop from the blast hit table and use that in the blastdbcmd (potentially recalculating them with + and - 200 bp )
Question: Hi everyone!!

I am new to text mining studies carried out in bioinformatics. I was learning how to use pubmed.mineR package for text mining using mulit-abstract txt file. I read a lot on RISmed package and easypubmed R packages but it seems they have some disadvantage when you have to use their output in pubmed.mineR.

pubmed.mineR uses a text file containing multiple abstracts in "abstract" format of Pubmed (one format abovet many other pubmed formats such as xml). easyPubmed have such function that can be exploited while RISmed doesn't it seems. On the other hand you can retrieve only 5000 abstracts at a time using easyPubmed while RISmed has no such upper limit.

I randomly chose a cancer type "oral cancer" and it had above 1 lakh PMIDs. RISmed successfully retrieved the abstracts(1.1 GB data) however it's output is incompatible for pubmed.mineR. On the other hand, though easypubmed had a compatible output, it has retrieval limit since it use PubMed API at the backend.

Is there a way using CLI to retrieve all ~1 lakh of abstracts in "abstract" format from pubmed since the website itself after June 2020 update has set a limit to 10,000 abstracts at a time. Here I attach a short code I used with easypubmed

library("easyPubMed")
search_topic <- 'oral cancer'
my_entrez_id <- get_pubmed_ids(search_topic)
my_entrez_id$Count
?fetch_pubmed_data
my_abstracts_txt <- fetch_pubmed_data(my_entrez_id, retmax = 142000, format = "abstract")
writeLines(my_abstracts_txt, con = "oral_cancer.txt")

Answer: You can use [Entrezdirect][1] (example output truncated). Be sure to sign-up for and use [NCBI_API_KEY.][2]

$ esearch -db pubmed -query "Oral cancer" | efetch -format abstract >> some_file

> 1. Am J Otolaryngol. 2020 Aug 16;41(6):102685. doi: 10.1016/j.amjoto.2020.102685.
> [Epub ahead of print]
>
> HPV vaccination practices and attitudes among primary care physicians since FDA
> approval to age 45.
>
> Petrusek J(1), Thorpe E(2), Britt CJ(2).
>
> Author information:
> (1)Department of Otolaryngology, Loyola University Medical Center, Maywood, IL
> 60153, United States of America. Electronic address: Jeffrey.Petrusek@lumc.edu.
> (2)Department of Otolaryngology, Loyola University Medical Center, Maywood, IL
> 60153, United States of America.
>
> PURPOSE: The aim of this study was to examine HPV vaccine administration
> practices since FDA approval to age 45 and assess knowledge regarding HPV and its
> association with oropharyngeal cancer.
> METHODS: A survey was distributed to 86 primary care physicians at Loyola
> University Medical Center. The survey contained 11 questions designed to capture
> HPV vaccination practices, knowledge of FDA approval, and barriers to
> vaccination.
> RESULTS: 46 (53%) physicians completed the survey and 45 responses were included.
> Among respondents who treat males ages 9-21 and females ages 9-26, the

[1]: http://bit.ly/entrez-direct
[2]: https://ncbiinsights.ncbi.nlm.nih.gov/2017/11/02/new-api-keys-for-the-e-utilities/
Question: I am new to bioinformatics and just for the sake of this question, suppose that I am looking for the full genome sequence of [*Escherichia coli*][1] and [*Saccharomyces cerevisiae*][2] in FASTA format. I thought that these would be fairly easy to find online.

For Escherichia coli I was only able to find some links to the ncbi website but I could not find anywhere in the website a link to download the genome in FASTA format. For Saccharomyces cerevisiae, I found [this FTP site][3] linked by ensembl with the sequence for each chromosome. I suppose that I could concatenate all the files in order to obtain the full genome in one file but isn't there already one file with the entire sequence? I saw the top level files but those only appear for masked and soft masked DNA, I think it would be doable to work with those if modern sequencers masked the DNA automatically, do they?

So, in summary, how do experienced people find full reference genomes for those two organisms and possibly other common species? I have only been using Google and the only page that appears to provide FASTA files for some species is Ensembl.

[1]: http://en.wikipedia.org/wiki/Escherichia_coli
[2]: http://en.wikipedia.org/wiki/Saccharomyces_cerevisiae
[3]: ftp://ftp.ensembl.org/pub/current_fasta/saccharomyces_cerevisiae/dna/

Answer: First place to check is Ensembl. They have separate pages for metazoa, fungi, etc. as well as the normal vertebrate page (which also contains S. cerevisiae). Each of the sites have a download page (click on "Downloads" at the top and then "Download data via FTP" on the right), where you can easily find fasta/GTF/etc. files.

BTW, a soft-masked genome is the genome. Soft-masking just means that predicted repeat regions are lower-case.

A second BTW is that you'll usually want the "toplevel" file.
Question: I have a single fasta file genome that contains 40Mbps spread throughout ~30,000 separate sequences (contigs). About half is expected to be repetitive DNA. I am looking for a tool to either:

1) cut repeats from the original file and paste to a new fasta file

2) delete repeat regions from file

or 3) mask repeat regions (replace all repetitive sequences with N)

The first option is ideal, but for any of the three choices I want to be as liberal as possible with the definition of "repetitive DNA". I want to avoid any potential repeat at all costs. Losing good data is better than keeping repeat data in this scenario.

Note that I don't want to reduce the number of times a sequence is repeated, but I want to delete or mask every instance of that repeat so that it is not found a single time in my genome file.

Any suggestions for tools that will perform any of these tasks? Thanks!

Answer: The R package Biostrings, from Bioconductor, will accomplish this. Here's a tutorial, including specific examples on masking. It's fast.

Edit: the example uses regions identified from RepeatMasker, if you haven't ID'd them yet.
Question: Hello!

I am very new to biopython and I am trying to accomplish what I think is a simple task: I would like to remove sequences from a protein alignment that do not contain a particular residue at a specified position. I would like to be able to input a protein alignment in fasta format and then output a new alignment where all the sequences that do not meet my criteria are removed

For example: My input protein alignment contains sequences that have a mixture of residues at position 137. I would like to output a new alignment that contains only sequences that have either an arginine or a valine at position 137.

Just a bit of additional clarification: I am sequencing an amplicon of a functional gene and generating protein sequence alignments using RDP's fungene pipeline. I want to further screen the alignment by eliminating any sequences that do not contain a selection of conserved residues at various positions.

Thank you very much for your time.

-J

Answer: Thank you so much for the answer!

It was just the advice I needed.

This is the python file I made to accomplish my goal (it's probably really ugly to anyone who has actual python experience)

from Bio import AlignIO

from Bio.Align import MultipleSeqAlignment
alignment = AlignIO.read("MyProt.fasta", "fasta") # my input alignment

goodseqs1 = MultipleSeqAlignment([]) # sets up an empty MSA that good sequences can be added to
goodseqs2 = MultipleSeqAlignment([]) # MSA that good seqs can be added to for another round of screening
badseqs = MultipleSeqAlignment([]) # MSA that seqs not meeting the criteria are added to
for sequence in alignment:
if sequence.seq[48] == "F":
goodseqs1.append(sequence) # adds all sequences with "F" at position 48 to goodseqs1 alignment
elif sequence.seq[48] == "Y":
goodseqs1.append(sequence) # adds all seqs with a "Y" at position 48 to goodseqs1 align
else:
badseqs.append(sequence) # puts all remaining seqs in badseqs

for sequence in goodseqs1: # additional round of screening for seqs that passed the first round
if sequence.seq[46] == "Q":
goodseqs2.append(sequence)
elif sequence.seq[46] == "R":
goodseqs2.append(sequence)
elif sequence.seq[46] == "G":
goodseqs2.append(sequence)
elif sequence.seq[46] == "I":
goodseqs2.append(sequence)
else:
badseqs.append(sequence)

AlignIO.write(goodseqs2, "SCREENED_SEQS.FASTA", "fasta") # writes a fasta alignment containing only passing seqs
AlignIO.write(badseqs, "Discard.fasta", "fasta") # writes a fasta alignment containing only failing seqs

print("Alignment length %i" % alignment.get_alignment_length()) # prints the length of the alignment
print("initial # of seqs", len(alignment)) # prints the # of seqs in the initial alignment
print("seqs passing first screen:", len(goodseqs1)) # prints the # of seqs passing the first screen
print("final # of seqs", len(goodseqs2)) # prints the # of seqs passing both screens
print("output files are: SCREENED_SEQS.FASTA and Discard.fasta") # prints the names of the output files
Question: I'm blasting+ (blastn+ or tblastn) an annotated sequence (a Genbank.xml sequence (nucleotide or protein) or an Uniprot.xml entry) against a DNA database.

Is there a standard tool to map the features of the query to the hit ? The best output would be a BED file.

Thanks.

Answer: I've answered my own question on my blog: http://plindenbaum.blogspot.fr/2013/07/mapping-annotations-of-query-sequence.html

$ java -jar dist/blastmapannots.jar I=P04514.xml B=blast_P04514_vs_AY065842.xml

AY065842 25 961 Non-structural_protein_3 943 + 25961 255,255,255 1 936 25
AY065842 34 469 RNA-binding 970 + 34 469 255,255,255 1 435 34
AY065842 472 640 Dimerization 947 + 472 640 255,255,255 1 168 472
AY065842 532 724 Interaction_with_ZC3H7B 917 + 532 724 255,255,255 1 192 532
AY065842 646 961 Interaction_with_EIF4G1 905 + 646 961 255,255,255 1 315 646
AY065842 520 733 coiled-coil_region 916 + 520 733 255,255,255 1 213 520
Question: Hi,

I have read that when one encounters a confounding variable, the preferred method to account for their effects would be to add this variable to the design formula but not remove them using limma's removebatcheffect function.
My questions on this topic are:
- should I remove batch effect using limma's removebatcheffect function if the confounding effect is large and inconsistent as suggested here (https://support.bioconductor.org/p/125386/#125387). For example, when one can see in a pca plot and heatmap that samples still cluster according to the presence or absence of a confounding variable even though it has been added to the design formula (eg when samples cluster by sex after it has been added to the design formula).

- should one use the transformed counts after removal of a confounding effect or the untransformed counts (and only add the confounding variable to the design formula) for heatmaps, pca plots and DE analysis?

Answer: Confounders should be modeled where that is possible, not removed. This particularly applies to DE analysis, or anything that relies on count-based statistics, like the Poisson or Negative Binomial distribution.

Where you are not modeling the data, but just presenting it - like in heapmaps or pca, then modeling the confounder isn't possible, and it should be removed.

Where the counfounder is known, you can directly use it in the design formula, or remove it using `limma::removeBatchEffects`. Where the confounder is unknown, then you must identify it using something like SVA.
Question: I like this dot plot (Figure a) for showing marker gene expression patterns for the clusters. Any ideas about how to make this plot with ggplot2? Some pseudocode would be enough.

![enter image description here][1]

[1]: https://media.springernature.com/full/springer-static/image/art%3A10.1038%2Fs41467-020-19737-2/MediaObjects/41467_2020_19737_Fig2_HTML.png?as=webp

Answer: Use `dittoDotPlot` from the [dittoSeq](https://bioconductor.org/packages/release/bioc/vignettes/dittoSeq/inst/doc/dittoSeq.html#561_dittoDotPlot) package. Much easier than from scratch. The resulting plot is still a ggplot2 object, so you can customize as you'd like. Though the function itself is also quite flexible. Works natively on SingleCellExperiment, Seurat, or SummarizedExperiment objects and can use gene expression and/or continuous metadata variables as desired.

So all you'd have to do is grab your top 5 markers for each cluster or what have you and feed them in along with your data object:

your_markers <- c("Gfap", "Vim", "Aqp4")
dittoDotPlot(sce_seurat_or_se_object, your_markers, group.by = "cluster_metadata_column") + coord_flip()
Question: We often run the Angilent 2100 after library construction. Our tutor thinks that the insert size should be controled at the range of 300bp to 400bp. I just don't know why the insert size has to e that. Hope someone to remove my confusion. Many thanks!

Answer: I think optimal insert size depends partly on the experimental design/objectives and sequencing chemistry used (e.g. MiSeq v3 PE 600 would be different than HiSeq PE 200), and partly due to clustering efficiencies. It doesn't _have_ to be 300-400 bp, but that might be the optimal range for the sequencing that you are going to perform.

Let's use this picture to help us:

![enter image description here][1]

Although it depends on your application, ideally you will have _some_ inner distance, avoiding overlap between Read 1 and Read 2. Let's say you didn't size select your library, and had inserts that ranged from 50 - 500 bp, and that you were sequencing your library on a HiSeq with 2 x 100 reads. After cluster generation, you would have a percentage of clusters that were < 200 bp in length, resulting in overlap of Read 1 and Read 2. Even worse, you'd have some clusters < 100 bp, meaning Read 1 and Read 2 would be identical and have some adapter read-through, "wasting" a certain amount of your sequencing power.

To compound the problem, smaller inserts tend to cluster more efficiently. From the Illumina Tech Note [Nextera Library Validation and Cluster Density Optimization][2]:

> It is important to consider library size when preparing samples for
> cluster generation. Because the clustering process preferentially
> amplifies shorter libraries in a mixture of fragments, large libraries
> tend to cluster less efficiently than small libraries.

So if you have a range of insert sizes, you expect the smaller inserts to "out-compete" the larger ones during cluster generation, making the above described issues more frequent. Instead, if you size-select your library to a certain range (say, 300-400 bp), you can ensure that your clusters will almost all have the same insert size, and that you will avoid overlap between Read 1 and Read 2 (if you use the right chemistry).

Why not select larger fragment sizes? As noted above, clustering efficiency decreases, and also larger fragments generate larger clusters, thus reducing the maximum cluster density of the flow cell.

[1]: http://www.frontiersin.org/files/Articles/77572/fgene-05-00005-HTML/image_m/fgene-05-00005-g001.jpg
[2]: https://www.illumina.com/documents/products/technotes/technote_nextera_library_validation.pdf
Question: Hi,

I wanted to query the gnomAD database without having to download and configure the entire data. Is there a way gnomAD provides an API for programmatic access? For example, for a given variant id, I want to retrieve the Total 'Allele Count':

INPUT: '1-978628-C-T' (or 'rs116586548')
OUTPUT: 180

i.e. if we query the website ([gnomAD][1]

[1]: https://gnomad.broadinstitute.org/variant/1-978628-C-T

I have thousands of variants that I want to get the Total Allele Count for, but would prefer not to download the entire dataset.

thanks!

Answer: use an API like https://myvariant.info/

wget -q -O - "http://myvariant.info/v1/variant/rs116586548" H "accept: application/json"
(..)
"an": {
"an": 30864,
"an_afr": 8690,
"an_amr": 838,
"an_asj": 302,
"an_eas": 1622,
"an_female": 13804,
"an_fin": 3490,
"an_male": 17060,
"an_nfe": 14940,
"an_oth": 982,
"an_popmax": 8690,
"an_raw": 30990
},

(...)
Question: Hi Biostars,

I have some prefiltered scRNAseq count data and want to extract counts for cells belonging to a specific cluster.

The following code retrieves the counts with gene symbols by default:

cluster01 <- subset(data_integrated,idents = c("Cluster 1"))
cluster01_counts <- cluster01@assays$RNA@counts

The output looks like this (column names truncated):

X1 AAAC.. AAAC.. AAAC..
A1BG 0 0 3
A1CF 0 0 0

Column X1 contains gene symbols.

Is there a way to extract ENSEMBL IDs directly instead of gene symbols?

Answer: This is what I would do in your situation.

# generate a table with ensemble and symbol for each gene in your data (good to keep it as a reference)
library(org.Hs.eg.db)
master_gene_table <- mapIds(org.Hs.eg.db, keys = cluster01_counts$X1, keytype = "SYMBOL", column="ENSEMBL")
master_gene_table <- as.data.frame(master_gene_table)

master_gene_table$ensembl <- cluster01_counts$ensembleID

Note that you need to make sure to remove duplicated genes if you want to transfer the ensemble names to `rownames(cluster01_counts)`