cannin/biostars_qa huggingface dataset:
Question: Hello, I have a file which gives me the position of the altered nucleotide or small insertion or deletion (reference allele and altered allele). Is there any tool which can be used to perform these changes on the reference nucleotide sequence and incorporate the ALT_ALLELE changes ?

CHROM	POS	REF_ALLELE	ALT_ALLELE	TRANSCRIPT_ID	GENE_ID	GENE_NAME	REGION	VARIANT_TYPE	REF_AMINO	ALT_AMINO	AMINO_POS	SIFT_SCORE	SIFT_MEDIAN	NUM_SEQS	dbSNP	SIFT_PREDICTION
SL3.0ch00	723860	A	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	NONSYNONYMOUS	W	G	52	0	4.32	1	novel	DELETERIOUS (*WARNING! Low confidence)
SL3.0ch00	723867	A	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	SYNONYMOUS	G	G	49	1	4.32	1	novel	TOLERATED
SL3.0ch00	723903	T	C	mRNA.Solyc00g005060.1.1	gene.Solyc00g005060.1	Solyc00g005060.1	CDS	SYNONYMOUS	G	G	37	1	4.32	1	novel	TOLERATED

Answer: Hello,

as https://www.biostars.org/u/40003/ says before [bcftools consensus][1] is the way you have to go. Therefor you need to convert your file to a valid `vcf` file.

$ echo "##fileformat=VCFv4.2" > input.vcf
$ echo "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO" >> input.vcf
$ tail -n+2 variants.txt|awk -v OFS="\t" '{print $1,$2,".",$3,$4,".",".","."}' >> input.vcf

Then compress it with `bgzip` and index by `tabix`.

$ bgzip -c input.vcf > input.vcf.gz
$ tabix input.vcf.gz

Now you are ready for `bcftools consensus`.

$ bcftools consensus -f genome.fa input.vcf.gz -o output.fa

fin swimmer

---
`bcftools consensus` may warn you about sequences not found in the `vcf`. You can safely ignore this. The warning happens if it finds sequence id's in the reference that are not in your vcf. If you want to get rid of those messages each contig of the reference needs to be defined in the vcf header. This can be done like this:

**1. Index the reference**

$ samtools faidx genome.fa

**2. Create headers for the `vcf`**

$ echo "##fileformat=VCFv4.2" > input.vcf
$ awk '{print "##contig="}' genome.fa.fai >> input.vcf
$ echo "#CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO" >> input.vcf

**3. Add variant information**

$ tail -n+2 variants.txt|awk -v OFS="\t" '{print $1,$2,".",$3,$4,".",".","."}' >> input.vcf

Continue with `bgzip`, `tabix` and `bcftools consensus` as above.

[1]: http://www.htslib.org/doc/bcftools.html#consensus
Question: If I run `fastq-dump --split-3` on a sra file, I get

```
file_1.fastq
file_2.fastq
file.fastq
```

My questions is how I handle `file.fastq`? Should I just ignore it?

Answer: `--split-3` will output 1,2, or 3 files: 1 file means the data is not paired. 2 files means paired data with no low quality reads or reads shorter than 20bp. 3 files means paired data, but asymmetric quality or trimming. In the case of 3 file output, most people ignore `.fastq`. This is a very old formatting option introduced for phase1 of 1000genomes. Before there were many analysis or trimming utilities and SRA submissions always contained all reads from sequencer. Back then nobody wanted to throw anything away. You might want to use `--split-files` instead. That will give only 2 files for paired-end data. Or not bother with text output and access the data directly using sra ngs apis.
Question: Hello,

I would like to fetch 1:1 orthologs for **GBRA2_HUMAN**, a cannonical GABAA receptor alpha-2 subunit sequence as defined by Uniprot (**Uniprot ID: P47869**), however it seems that this sequence is not available in the OMA database. Could you please let me know if there is an alternative way to obtain the orthologs for this gene and what might be the reason for the missing entry?

I tried using the **OMAdb** R package:

getHOG(id ='GBRA2_HUMAN', members = TRUE)

>THE OMA REST API request failed:https://omabrowser.org/api/hog/GBRA2_HUMAN/members/
Here's the original error message:
Not Found (HTTP 404).

If I use `mapSequence(sequence)` using the cannonical sequence from Uniprot, I get three matched targets, neither of which is the cannonical human GABAA alpha2 subunit:
>[1] "A0A2R8ZZ28" "H2QPE5" "G3QVC6"

On the other hand, other human GABAA subunits seem to be present in the database (e.g. **GBRA1_HUMAN**, **GBRA3_HUMAN**).

The **GBRA2_HUMAN** sequence I used for searches is the following:

sequence <- 'MKTKLNIYNMQFLLFVFLVWDPARLVLANIQEDEAKNNITIFTRILDRLLDGYDNRLRPG
LGDSITEVFTNIYVTSFGPVSDTDMEYTIDVFFRQKWKDERLKFKGPMNILRLNNLMASK
IWTPDTFFHNGKKSVAHNMTMPNKLLRIQDDGTLLYTMRLTVQAECPMHLEDFPMDAHSC
PLKFGSYAYTTSEVTYIWTYNASDSVQVAPDGSRLNQYDLLGQSIGKETIKSSTGEYTVM
TAHFHLKRKIGYFVIQTYLPCIMTVILSQVSFWLNRESVPARTVFGVTTVLTMTTLSISA
RNSLPKVAYATAMDWFIAVCYAFVFSALIEFATVNYFTKRGWAWDGKSVVNDKKKEKASV
MIQNNAYAVAVANYAPNLSKDPVLSTISKSATTPEPNKKPENKPAEAKKTFNSVSKIDRM
SRIVFPVLFGTFNLVYWATYLNREPVLGVSP'

Many thanks.

Answer: The reason for this is that OMA is quite stringent when mapping IDs. For human, OMA uses the genome provided by Ensembl. The protein sequence we consider for the gene you are interested in (human gene "GABRA2") is a different isoform in OMA vs UniProt. The one we use in OMA (https://omabrowser.org/oma/info/ENSP00000427603) which is 511 AA long and maps to a different UniProt entry (https://uniprot.org/uniprot/E9PBQ7).

So how can you find GBRA2_HUMAN in OMA?

1) One way would be to use the approximate search functionality, e.g. using the REST API function sequence list (setting the search parameter to "approximate").
> https://omabrowser.org/api/docs#sequence-list

For instance, using the library [OmaDB][1] in R:
`mapSequence(sequence,search = "approximate")`

2) Another solution would be to map via the gene name GABRA2. You can use the REST API function xref list:
>https://omabrowser.org/api/docs#xref

For instance, using the library [OmaDB][1] in R:
`searchProtein('GABRA2')`

This returns two human isoforms in OMA, the first of which is ENSP00000427603 (and which was used to infer orthologs).

[1]: https://bioconductor.org/packages/release/bioc/html/OmaDB.html
Question: Hi there,

I know there are different ways to extract some "wanted" sequences from multifasta file. In general, having one or more IDs, you can retrieve the sequence from multifasta file using bioperl script or another kind of software as ones pointed out here (https://www.biostars.org/p/2822/). In my case I know that I want the sequence number 1, 5 and 10 from a multifasta file. Any idea how can I achieve this using perl/bioperl preferably?

Thanks in advance.

P.D: I'm a beginner in bioinformatics word, sorry if it is a very basic question.

Answer: Here is a BioPerl example since you asked:

https://gist.github.com/sestaton/651a34ae8dbb4d4c945c

I put this in a script so it is easier to read, but in practice, I would just do this at the command line with a one-liner (that's how I wrote it):

perl -MBio::SeqIO -M'List::MoreUtils qw(any)' -E '@ids = (1, 5, 10); $num = 0; $seqio = Bio::SeqIO->new(-fh => \*STDIN); while ($seq = $seqio->next_seq) { $num++; say join "\n", ">".$seq->id, $seq->seq if any { $_ == $num } @ids; }' < seqs.fas

That will save you from opening an editor. It's not as readable but it works the same. :)
Question: **Hi Everyone**

I am trying to use Michigan Impute Server. I use Checkvcf first to avoid failure in the server.

python checkVCF.py -r checkVCF/hs37d5.fa -o test chr3.vcf

**I got some duplicates and inconsistent ref.**

> checkVCF.py -- check validity of VCF file for meta-analysis version
> 1.4 (20140115) contact zhanxw@umich.edu or dajiang@umich.edu for problems. Python version is [ 2.7.5.final.0 ] Begin checking vcfFile
> [ chr3.vcf ] Duplicated site [ 3:14187449 ] Duplicated site [
> 3:21307401 ] Duplicated site [ 3:38608045 ] Duplicated site [
> 3:39146429 ] Duplicated site [ 3:41912651 ] [ 10000 ] lines processed
> Duplicated site [ 3:48618728 ] Duplicated site [ 3:79399575 ]
> Duplicated site [ 3:95176677 ] Duplicated site [ 3:96472739 ]
> Duplicated site [ 3:99067458 ] [ 20000 ] lines processed Duplicated
> site [ 3:113876275 ] Duplicated site [ 3:120522716 ] Duplicated site [
> 3:121633904 ] Duplicated site [ 3:128622922 ] [ 30000 ] lines
> processed Duplicated site [ 3:171926373 ] Duplicated site [
> 3:183371250 ]
> --------------- REPORT --------------- Total [ 37146 ] lines processed Examine [ 33 ] VCF header lines, [ 37113 ] variant sites, [
> 378 ] samples [ 16 ] duplicated sites [ 0 ] NonSNP site are outputted
> to [ test.check.nonSnp ] [ 6995 ] Inconsistent reference sites are
> outputted to [ test.check.ref ] [ 0 ] Variant sites with invalid
> genotypes are outputted to [ test.check.geno ] [ 0 ] Alternative
> allele frequency > 0.5 sites are outputted to [ test.check.af ] [ 0 ]
> Monomorphic sites are outputted to [ test.check.mono ]
> --------------- ACTION ITEM ---------------
> * Remove duplicated sites and rerun checkVCF.py
> * Read test.check.ref, for autosomal sites, make sure the you are using the forward strand
> * Upload these files to the ftp server (so we can double check): test.check.log test.check.dup test.check.noSnp test.check.ref
> test.check.geno test.check.af test.check.mono

How can I remove this **duplicate sites and inconsistent reference sites**?

I tried this but it seems it excludes duplicate variants not sites:

plink --bfile snps_filtered --list-duplicate-vars ids-only suppress-first
plink --bfile snps_filtered --exclude plink.dupvar --make-bed --out snps.DuplicatesRemoved
plink --bfile snps_filtered --recode vcf --snps-only just-acgt --out snps.final

A link to where is this in plink will be okay too.

**Thanks**

Answer: I note that you are comparing to *hs37d5.fa*, but to which genome was the original sample aligned? It would likely be recorded in your VCF header.
Question: I re[cently asked about help finding promoters in particular][1], and got a helpful pointer for that... however, it wasn't exactly what I wanted I found out.

I then found something much closer to what I wanted: the [Integrated Genome Viewer has a "Motif Finder"][2] tool, and you can even enter in ambiguous nucleotides to it. The only catch (for me at least) is... is there a command line tool like it? Do I have to interact with IGV through the GUI? It seems there is some command line support for IGV, [but it seems like it is geared more for pre-processing][3]. Sorry if I missed something though, would appreciate any pointers!

[1]: http://%20https://www.biostars.org/p/484669/#484737
[2]: https://software.broadinstitute.org/software/igv/motif_finder
[3]: https://software.broadinstitute.org/software/igv/igvtools_commandline

Answer: If any useful, I wrote [fastaRegexFinder](https://github.com/dariober/bioinformatics-cafe/tree/master/fastaRegexFinder) for this.
Question: This question was inspired by Peter Cock 's tweet

https://twitter.com/pjacock/status/118012750105546752

Most softwares (C, java...) use a int32 (unsigned or not) to store the length of the chromosome. It isn't enough when the length of a chromosome is greater than INT_MAX or UINT_MAX

# define INT_MAX 2147483647
# define UINT_MAX 4294967295U

So my question is:

Is there any resource where one can find the length of the chromosomes. Something like BioNumbers.

What's the length of the longest chromosome ?

Answer: [This old poster][1] gives highest genome size of 130 billion base pairs for lungfish

Biochemical measurements [suggest that certain Amoeba might have larger genomes][2]

Also see references for [this article][3] and [this one][4].

(note that these are genome sizes not chromosome sizes but large genome size and low chromosome count might mean chromosomes are large assuming most of genome isn't mitochondrial DNA or something)

[1]: http://www.sciencemag.org/site/products/posters/SequencingPoster.pdf
[2]: http://www.genomesize.com/statistics.php
[3]: http://en.wikipedia.org/wiki/Genome_size#Conversion_from_picograms_.28pg.29_to_base_pairs_.28bp.29
[4]: http://en.wikipedia.org/wiki/Genome#Genome_size
Question: Hi,

I just had a quick question that I started to confuse myself with. When I am about to index my reference sequence for some RNA-Seq data and I got confused if I use the cDNA or CDS or transcript fasta file.

In the past I used a transcript fasta but I can't seem to find one for the species an am working with from ensembl (http://useast.ensembl.org/info/data/ftp/index.html/)

Maybe I can just use the whole genome for the reference? I am using Salmon to map the reads.

Thanks in advance

Answer: Here is an [illustration of differences][1] between CDS/transcripts.

You can find cDNA's and CDS's for *Felis catus* [here][2].

Ideally you should align your data to the genome and then use an annotation file to do counting of reads.

Added: If you are using `salmon` then you should use transcripts.

[1]: https://biology.stackexchange.com/questions/11377/difference-between-cds-and-cdna
[2]: ftp://ftp.ensembl.org/pub/release-92/fasta/felis_catus/
Question: I am working on the Differential Expression Analysis on RNA-Seq Data between tumor sample and normal sample, while my lecturer suggest me better use the raw expression data instead of using FPKM as expression unit , anyone could tell me why ? or what is the advantage for this ?

Answer: In a nutshell, because you'll be able to apply statistical models to your values. When you normalize to FPKM you lose your ability to do so and end up with less robust statistical tests
Question: Im trying to make a 3D PCA but i m not sure if its really working or not i cannot label my sample which i intend to here is my code anyone can help or any suggestion how can i label my samples in my final plot.

d <- read.csv("~/analysis/gene_fpkm/PRIMARY_AUX_ATPASE.csv",header = TRUE)
names(d)
model <- prcomp(t(d[,-1]), scale=TRUE)
plotPCA <- function(x, nGroup) {
n <- ncol(x)
}
fit <- hclust(dist(x), method="complete")
groups <- cutree(fit, k=nGroup)
if(n == 3) { # 3d plot
plot3d(x, col=groups, type="s", size=1, axes=F)
axes3d(edges=c("x--", "y--", "z"), lwd=3, axes.len=2, labels = TRUE)
mtext3d(text)
grid3d("x")
grid3d("y")
grid3d("z")
} else { # 2d plot
maxes <- apply(abs(x), 2, max)
rangeX <- c(-maxes[1], maxes[1])
rangeY <- c(-maxes[2], maxes[2])
plot(x, col=groups, pch=19, xlab=colnames(x)[1], ylab=colnames(x)[2], xlim=rangeX, ylim=rangeY)
lines(c(0,0), rangeX*2)
lines(rangeY*2, c(0,0))
}
}
plotPCA(model$x[,1:3], 5)

this is the code im running i cannot label my sample in which comes in the final 3D pca .any suggestion or help would be highly appreciated

Answer: This is a question more suited to an R programming forum. Anyway, I suggest you look into the [car package][1] and its scatter3d() function (hint: look for the labels parameter). Otherwise, with rgl, you need to use the text3d() function.
The statistical tools for high-throughput data analysis (STHDA) web site has [tutorials on 3D graphics in R][2] (among other useful tutorials).

[1]: https://CRAN.R-project.org/package=car
[2]: http://www.sthda.com/english/wiki/3d-graphics
Question: Hi everyone,

I’m analyzing a RNAseq dataset and I’ve come across something that is giving me some headaches.

First some background: I have eight cyanobacterial samples from two different strains (four vs four) that show a (rather small) phenotypic difference under the same experimental conditions. My goal is to identify differentially expressed genes that can potentially explain this phenotypic difference. Though I wasn’t expecting to find a lot of DGEs since the strains are very similar.

Now the data; After aligning and counting, I found the total number of counts differs between samples (figure 1A), as it does too the ribosmal RNA depletion efficiency (figure 1B). All normal up to here, the rRNA depletion kits are not optimized for cyanobacteria and it is not the first time I see this discrepancy in efficiency between samples. A couple more observations, the percentage of rRNA correlates with the total number of counts (figure 1C) and after removing the rRNA counts, all samples have similar number of counts (figure 1D).

![figure_1][1]

I removed the rRNA counts and used DESeq2 to perform variance stabilization on the counts and perform PCA to see how the samples were clustering (figure 2). PC1 explains ~ 50% of the variance and, at first glance, it Iooked like it was separating the strains with the exception of S4. An outlier I thought.

![figure_2][2]

But then I looked again at the percentage of rRNA in each sample and realized the samples’ clustering in PC1 could also be explained by the rRNA percentage, I plotted this and the correlation is there (figure 3).

![figure_3][3]

So I’m wondering a few things.

1. Does it make any sense that most of the variation in the data is explained by the percentage of rRNA?
2. Or is this just an artifact?
3. If it is not, how can this be explained?
4. Does it make the data unusable for my original question?

Thank you in advance!

Note: in case it matters, the ribosomal depletion was performed with MICROBExpress kit.

[1]: /media/images/6fc7f854-0d9b-4b80-91c2-cf436c5d
[2]: /media/images/f3b8b3e4-ba7a-43dc-82a9-8eab3bdb
[3]: /media/images/b22c0ed1-3261-4c7c-beab-dcc07bbe

Answer: As a practical note: If you worry that rRNA proportion impacts the estimation of other genes, you can (and likely should) include it as a covariate in DESeq2 with a design `~ strain + rrna_proportion`.

As a technical note: it seems like you removed rRNA counts *prior* to normalization -- meaning that the library sizes associated with the counts are possibly not correct. The easiest way to adapt your current workflow would be to run `estimateSizeFactors` on the original counts; then remove ribosomes, and then use `sizeFactors(de.obj)` to pass in the originally-estimated size factors. This could result in a consistent over-estimation of expression of non-ribosomal genes; resulting in the observed PC1.
Question: I have a small dataset containing five genes and their expression at 20 different time points, I want to construct the transcriptional regulatory network based on the data. Are there any tools or tutorials available to do this? I would much prefer tutorials.

Answer: I was finally successful to do this using both of the following approaches, I found the first one to be easier to work with:

1. Äijö T., Lähdesmäki H. (2009) Learning gene regulatory networks from
gene expression measurements using non-parametric molecular
kinetics. Bioinformatics, 25, 2937–2944.
http://www.cs.tut.fi/~aijo2/gp/

2. Bonneau R, Reiss DJ, Shannon P, et al. The Inferelator: an algorithm
for learning parsimonious regulatory networks from systems-biology
data sets de novo. Genome Biology. 2006;7(5):R36.
doi:10.1186/gb-2006-7-5-r36.
http://baliga.systemsbiology.net/drupal/content/inferelator/
Question: Maybe I'm missing something, but why there are exons that the start and end are the same?!

Examples:

1	havana	exon	2768091	2768091	.	+	.	gene_id "ENSG00000285945"; gene_version "1"; transcript_id "ENST00000648684"; transcript_version "1"; exon_number "1"; gene_name "AC242022.2"; gene_source "havana"; gene_biotype "lncRNA"; transcript_name "AC242022.2-201"; transcript_source "havana"; transcript_biotype "lncRNA"; exon_id "ENSE00003831437"; exon_version "1"; tag "basic";

17	ensembl	exon	41624191	41624191	.	-	.	gene_id "ENSG00000128422"; gene_version "17"; transcript_id "ENST00000540235"; transcript_version "5"; exon_number "2"; gene_name "KRT17"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "KRT17-205"; transcript_source "ensembl"; transcript_biotype "protein_coding"; exon_id "ENSE00003733771"; exon_version "1"; tag "basic"; transcript_support_level "5";

Thanks

Answer: Hmmm seems like an off-by-one error, as these exons are actually a single base long. But GTF is one-based so the coordinates indicate, indeed, a single nucleotide: https://www.biostars.org/p/84686/
Question: Hi there,

We are using Nanopore for some transcriptomics analysis and I am building a QC pipeline. There is a default statement that it could be good to remove reads below 500 bp in length but I could not find arguments for or against this in the literature.

So my question is, assuming that the quality of the read is good, is there an argument why one should remove shorter reads below a certain threshold ( except of course like super short stuff like 10 or so) ?

Thanks!

Answer: As always, one needs to know which downstream analysis you have in mind.

It is recommended, yes, but it doesn't make much sense for transcriptomic data. Many transcripts are [well below 500 nt][1] in length, so you tend to lose information. I lately noticed this in `BUSCO`-analyses. My suggestion would be to try length filtering, and compare it to other measures (while watching out for transcriptome completeness). Filtering for quality with [nanofilt][2], and selection of full-length transcripts using something like [pychopper][3] would be more appropriate in my opinion.

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7905317/
[2]: https://github.com/wdecoster/nanofilt
[3]: https://github.com/nanoporetech/pychopper
Question: Hello,

I can go through the Seurat all the code from here: https://satijalab.org/seurat/archive/v3.0/visualization_vignette.html.

My question is: Is that possible that I can get the frequency of one gene in one cluster?

Thank you in advance for your great help!

Best Regards,

Yue

Answer: `DotPlot` computes the fraction of cells expressing a gene in each group (metadata column), so you can cheat a little and extract the results from it [[source](https://github.com/satijalab/seurat/issues/1888)].

library("Seurat")

features <- c("PECAM1", "CD14")
group_column <- "clusters"

perc_exp <- DotPlot(seu, features=features, group.by=group_column)$data[, c("features.plot", "id", "pct.exp")]

The returned data will look something like this:

> perc_exp
features.plot id pct.exp
1 PECAM1 Endothelial 26.55087
2 CD14 Endothelial 37.21239
3 PECAM1 T Cells 57.28534
4 CD14 T Cells 90.82078

You can then plot this easily using ggplot.

library("ggplot")

ggplot(perc_exp, aes(x=id, y=pct.exp, fill=id)) +
geom_col() +
facet_wrap(~features.plot)

![enter image description here][1]

----------

Example data used.

perc_exp <- structure(list(features.plot = c("PECAM1", "CD14", "PECAM1",
"CD14"), id = c("Endothelial", "Endothelial", "T Cells", "T Cells"
), pct.exp = c(26.55086631421, 37.212389963679, 57.2853363351896,
90.8207789994776)), class = "data.frame", row.names = c(NA, -4L
))

[1]: /media/images/9474f865-212b-4bbf-beb4-00fa57a5
Question: Hi,

I have many single sample vcf files containing structural variants. I'm trying to merge all the vcfs to one multi-sample vcf, however as the breakpoints of structural variants are only estimates, the position of the same structural variant in different samples would be slightly different. Is there a way to merge individual vcfs containing structural variants into one multi sample vcf? taking account for example a percentage of overlap to call structural variants in different samples as the same variant.

Thanks

Answer: For that purpose I use [SURVIVOR merge][1], which has a distance setting specifying the maximal difference between pairs of breakpoints (start1 vs start2 and end1 vs end2). As such it doesn't use a percentage of overlap, but it is also not requiring exact breakpoints.

About the percentage overlap the author of SURVIVOR said that he doesn't think it is correct that large SVs easier reach that minimal percentage, and therefore he just uses a distance. I believe he uses a distance of 1kb, but I commonly use 500bp.

[1]: https://github.com/fritzsedlazeck/SURVIVOR