cannin/biostars_qa huggingface dataset:
Question: Hi,

I was able to plot the PCA between two different variables using function autoplot() function from ggfortify library. I have a question about how to add the black borders around the circles and triangles as shown in the plot below. Please let me know if there any specific function for adding?

Thank you,

Toufiq

## PCA on two variables##
pca_v1 <- prcomp(Node, scale. = TRUE)
autoplot(pca_v1,
data = Neg_Dct, colour = 'Category', shape = 'Response' , size = 3.5) +
theme_classic()

![](https://i.imgur.com/RPst3AR.png)

Answer: We can achieve the same using basic *ggplot*, see example:

library(ggplot2)

#example data
df1 <- data.frame(
x = 1:10,
y = 1:10,
grp = rep(LETTERS[1:2], 5)
)

ggplot(df1, aes(x, y, col = "black", fill = grp, shape = grp)) +
geom_point() +
scale_color_identity() +
scale_fill_manual(values = c("red", "green")) +
scale_shape_manual(values = c(21, 24)) +
theme_classic()

We need to use shapes that can have colour and fill separately: http://www.sthda.com/english/wiki/r-plot-pch-symbols-the-different-point-shapes-available-in-r
Question: *Hello everyone*!
My name is Kseniia and I am new here!

I did encounter some problems while analyzing american gut data. Here is the code for it:

library(devtools)
devtools::install_github("zdk123/SpiecEasi")
library(SpiecEasi)
SpiecEasi.matrix <- spiec.easi(amgut1.filt,
method = 'glasso',
lambda.min.ratio = 1e-2,
nlambda = 20,
icov.select.params=list(rep.num=50))

The next line suppose to **add OTU names to rows and columns**

rownames(SpiecEasi.matrix$refit) <- colnames(amgut1.filt)

Then this error appeared:
**Error in `rownames<-`(`*tmp*`, value = c("OTU_1", "OTU_2", "OTU_3", "OTU_4", : attempt to set 'rownames' on an object with no dimensions**

I understood that it appears because $refit is a list. But I don't know how to fix it. If anyone had the same error, I would greatly appreciate your help with its solution.

Thanks in advance!
Kseniia

Answer: Try this:

rownames(SpiecEasi.matrix$refit$stars) <- colnames(amgut1.filt)
Question: I have multiple single-sample VCF files, which I want to merge into a single multi-sample VCF file. When using `bcftools merge` I am getting duplicate records.

$ bcftools merge ms01e_phased.vcf.gz ms02g_phased.vcf.gz ms03g_phased.vcf.gz ms04h_phased.vcf.gz MA605_phased.vcf.gz MA611_phased.vcf.gz -O v -o RBphased_variants.SixSamples.Final.vcf
# duplicate records at the same lines from the file "RBphased_variants.SixSamples.Final.vcf"
2	14691373	.	A	.	1153.31	PASS	BaseQRankSum=2.02;ClippingRankSum=0;ExcessHet=3.0103;FS=1.098;InbreedingCoeff=-0.0861;MQ=58.74;MQRankSum=-2.459;QD=19.22;ReadPosRankSum=-0.466;SOR=0.96;DP=68;AN=8	GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC:PM	0/0:4:4:0:0:0/0:.:.:0/0:.:.	0/0:7:7:18:0:0/0:.:.:0/0:.:.	0/0:4:4:12:0:0/0:.:.:0/0:.:.	0/0:2:2:3:0:0/0:.:.:0/0:.:.	./.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.
2	14691373	.	A	AAG	1153.31	PASS	BaseQRankSum=2.02;ClippingRankSum=0;ExcessHet=3.0103;FS=1.098;InbreedingCoeff=-0.0861;MQ=58.74;MQRankSum=-2.459;QD=19.22;ReadPosRankSum=-0.466;SOR=0.96;set=InDels;DP=676;AF=0.042;AN=4;AC=0	GT:AD:DP:GQ:PGT:PID:PL:PG:PB:PI:PW:PC:PM	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	0/0:12,0:12:9:.:.:0,9,135:0/0:.:.:0/0:.:.	0/0:22,0:22:12:.:.:0,12,180:0/0:.:.:0/0:.:.
2	14691374	.	A	.	1320.25	PASS	BaseQRankSum=-1.049;ClippingRankSum=0;ExcessHet=0.2929;FS=0;InbreedingCoeff=0.4006;MQ=55.35;MQRankSum=0;QD=33.01;ReadPosRankSum=-0.671;SOR=0.892;DP=44;AN=2	GT:AD:DP:GQ:PL:PG:PB:PI:PW:PC:PM	0/0:4:4:0:0:0/0:.:.:0/0:.:.	./.:7:7:.:0:./.:.:.:./.:.:.	./.:0:0:.:0:./.:.:.:./.:.:.	./.:0:0:.:0:./.:.:.:./.:.:.	./.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.
2	14691374	.	A	G	1320.25	PASS	BaseQRankSum=-1.049;ClippingRankSum=0;ExcessHet=0.2929;FS=0;InbreedingCoeff=0.4006;MQ=55.35;MQRankSum=0;QD=33.01;ReadPosRankSum=-0.671;SOR=0.892;set=HignConfSNPs;DP=710;AF=0.115;MLEAC=3;MLEAF=0.115;AN=4;AC=0	GT:AD:DP:GQ:PGT:PID:PL:PG:PB:PI:PW:PC:PM	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	./.:.:.:.:.:.:.:.:.:.:.:.:.	0/0:12,0:12:9:.:.:0,9,135:0/0:.:.:0/0:.:.	0/0:22,0:22:12:.:.:0,12,180:0/0:.:.:0/0:.:.

I raised this issue in `bcftools` thinking if it was a bug https://github.com/samtools/bcftools/issues/754 . But, is there any other solution to the problem.

Answer: Hey kirannbishwa01,

I would not call these duplicate records as they are different calls but at the same position.

Just adding `--merge all` to your command should help to solve it

So:

bcftools merge --merge all ms01e_phased.vcf.gz ms02g_phased.vcf.gz ms03g_phased.vcf.gz ms04h_phased.vcf.gz MA605_phased.vcf.gz MA611_phased.vcf.gz -O v > RBphased_variants.SixSamples.Final.vcf

Kevin
Question: I have a ligand that exists in its own PDB (*reference PDB*). I have then generated another (structurally similar, but not exact) ligand from SMILES, and written that to a PDB.

I am trying to find out how to superimpose the generated ligand onto the reference PDB, but I am having trouble finding a simple way to do this. This is a step that is part of a pipeline, so I need to be able to do it via command line, however, the tool I use can be flexible. I have been looking at `rdkit`, but am having trouble with that (the number of atoms don't match since the ligand compounds are similar, and not identical). It looks like PyMol can do this, but I don't know how to do it from the command line.

**Does anyone have any idea how to superimpose two PDB ligands from two PDB files?**

**UPDATE:**

The ligand generated from SMILES ligand energy minimized so it **does** have proper 3d coordinates.

Answer: I think I hav found an answer using OpenBabel.

obfit "Oc1ccc(cc1O)C1Nc2cccc3cccc(N1)c23" fixed.pdb moving.pdb > new_aligned.pdb

Where the smiles pattern matches my moving.pdb (which is where my generated conformers came from).
Question: Dear All,

I am having trouble with P value conversion for my data. I tried following [this post][1] but could not understand how I can get the P value out of LogP value or -log10 (p-value). Please share your experience I can get P value out of the data in either R or excel.

```
rsid CHR LOGPVALUE
rs3243810 1 3.835343448
rs3287090 1 3.835343448
rs3283061 1 2.225202326
rs3393481 1 2.225202326
```

Regards
Devender

[1]: https://www.biostars.org/p/358046/

Answer: Hi,

To obtain the *p*-value from the **`-log10(p-value)`** you can try (assuming that you are working from **R**, and the example data that you presented is called `data`):

data$pvalue <- 10**-(data$LOGPVALUE)

You can apply the formula on excel too: https://stackoverflow.com/questions/53913175/how-to-do-an-inverse-log-transformation-in-r

I hope this helps,

António
Question: Hi all,

I have DNA sequencing data where human tumors have mouse reads as per the experiment protocol(similar to PDX).

I have calculated the read counts of these mouse genes. using HTSEQ-count,after this I want to select only those mouse genes which have no homology to the human genome and only mouse specific.

Is there any tool to do this besides BLAST?Let me know.

Thanks,

Ron

Answer: Find all the orthologous genes between the two species (one database of many here, one how-to here), the mouse genes not present in this list are the genes you want.
Question: I would like to get non-variant sites using bcftools but am getting an odd error about the DP field.

bcftools mpileup -f genomic.fa file.bam | bcftools call --gvcf 5 -m > gvcf.vcf

--gvcf output mode requires FORMAT/DP tag, which is not present in the input header

Yet, the DP field is clearly in the output of mpileup:

NW_006532014.1 14714 . A <*> 0 . DP=52;I16=26,17,0,0,1926,98680,0,0,2580,154800,0,0,750,16322,0,0;QS=1,0;MQSB=1;MQ0F=0 PL 0,129,255
NW_006532014.1 14715 . T <*> 0 . DP=53;I16=27,18,0,0,1959,97077,0,0,2700,162000,0,0,794,17716,0,0;QS=1,0;MQSB=1;MQ0F=0 PL 0,135,255
NW_006532014.1 14716 . A <*> 0 . DP=50;I16=25,19,0,0,1905,91787,0,0,2640,158400,0,0,790,17884,0,0;QS=1,0;MQSB=1;MQ0F=0 PL 0,132,255
NW_006532014.1 14717 . T <*> 0 . DP=50;I16=25,18,0,0,1841,88215,0,0,2580,154800,0,0,766,17490,0,0;QS=1,0;MQSB=1;MQ0F=0 PL 0,129,255

Any idea what is the problem? When getting variant sites with `-vm` there is no issue. What I want is the DP4 field on a gvcf file, are the relevant numbers the first 4 after `I16=` ?

Thanks!

Answer:

bcftools mpileup -f genomic.fa file.bam -a "AD,ADF,ADR,DP,SP,INFO/AD,INFO/ADF,INFO/ADR" | bcftools call --gvcf 5 -m > gvcf.vc

Answer: https://samtools.github.io/bcftools/bcftools.html#mpileup

try to specify `--gvcf 5` and `--annotate` in **mpileup**

> annotate: Comma-separated list of FORMAT and INFO tags to output. (case-insensitive, the "FORMAT/" prefix is optional, and use "?" to list available annotations on the command line) [null]:
Question: Given a sam file which was produced by using bwa to align a fastq of reads to a reference genome, how can I enumerate the number of reads which were assigned to each chromosome of the reference genome?

Answer: Here's how I count the number per sequence of primary alignments with proper read pairs:

https://gist.github.com/slowkow/52a8911e18a873a51dcb
Question: Hi there,

there might be a rather simple solution to this but I already checked NCBI for a script (or something like this) but I wasn't lucky. So here is my problem...

I got a list of several species names, say 'Gardnerella vaginalis' is one of them, from that name I need to conclude what the genus is this species falls in. In most cases this is fairly easy since naming convention says that the first part is genus and second part specifies the species. However, since there might be one or the other synonym for a species (especially in the fungal world) I need some kind of lookup into the ncbi taxonomy which gives me the respective genus name for a given species name. (and in the best case this sould be script-able)

Any ideas are highly appreciated!

Thanks,

phil

Answer: With Entrez Direct:

```
esearch -db taxonomy -query "Gardnerella vaginalis" | efetch -format docsum | xtract -element Genus
Gardnerella
```

Alternatively, e.g.:

```
esearch -db taxonomy -query "Gardnerella vaginalis" | efetch -format xml | xtract -element Lineage
cellular organisms; Bacteria; Actinobacteria; Actinobacteria; Actinobacteridae; Bifidobacteriales; Bifidobacteriaceae; Gardnerella
```
Question: When I use samtools sort bam file, do i need to use samtool index the sorted bam file?
what kind situation i need to use samtools index?
Thanks!

Answer: > when i get the mapped sam file, i use samtools view to transfer this
> sam file to bam file, then i use samtools sort this bam file

Based on your explanation you do, the following, or something similar:

samtools view myfile.sam -o myfile.bam
samtools sort myfile.bam -o myfile_sorted.bam

That's not wrong, but it's also not necessary. Note that you can do the following in one go:

samtools sort myfile.sam -o myfile_sorted.bam

Finally, often you can also have your aligner write directly to samtools sort:

bwa mem genome.fa reads.fastq | samtools sort -o myfile_sorted.bam

As such you can avoid quite some wasted time and intermediate files.

> Depending on your reply, if i do not index the bam file it will not
> effect my later analysis( first condition), do i understand right?

A bam is a binary (compressed) format. For some applications, e.g. visualization in IGV or variant calling the tool needs 'random access' to the bam file, it needs to be able to find reads without iterating through the entire file. That's what the index is for.

We don't know what your "later analysis" is (please try to be as informative as possible when asking questions), but either it will work (no index required) or it will fail (and tell you it needs an index). Index files are generally small and created quickly, so there is no need to worry about them. Just index your bam.

> How to define " coordinate sorted"?

A bam is coordinate sorted if the reads are sorted by coordinates. So reads from the beginning of the first chromosome are first in the file. That's what `samtools sort` is for. The opposite would be 'sorted by name' in which reads are sorted by their read ID. Aligners also generate unsorted bams.
Question: I want to keep track of how many reads are overlapping each nucleotide in my genome, so that I can quickly get the number of overlapping reads with dict[chromosome][strand][nucleotide].

This is possible with intervaltrees from bx-python, but afair, that module did not work with multiprocessing, and it might not be the fastest for this specific task.

Are there other modules/packages that do this?

Answer: With python you have a few options depending on your expertise. The numbering does not reflect my perferences. The ease and speed of each option I listed is dependent on the size of your input and the number of regions to parse. If you have a lot of regions I would recommend using samtools.

1. [pybedtools][1] is a great wrapper of Bedtools in python. It can get slow depending on how many times you call the functions and how large your data are. However if you're familiar with intersection positions in Bedtools this should be a breeze!
2. Cython. You can use C types and np.ndarrays to iterate much faster in Cython
3. Use samtools: For your question
> how many reads are overlapping each nucleotide

`samtools depth` is perfect for this!

4. pysam: same thing a 3. But now you are using a htslib API in python.

[1]: https://pythonhosted.org/pybedtools/
Question: Hello,

I have performed a motif enrichment analysis with MEME-ChIP.

I am getting a beautiful motif with an e-value of 1.3e-357. And luckily I get a known motif (show more section below the logo) with an e-value of 1.1e-264.

When I click on "MEME" under the "Discovery/​Enrichment Program" section, I can see that 306 sites contributed to build the motif.

Knowing that I submitted 2359 sequences, I would tend to say that this motif is not significant.

Do you agree on that?

edit: The known motif which was found with centrimo, indicates 1658 in "region matches" section. does it mean that 1658 of my sequences have the motif?

Thanks a lot

Answer: Reply from MEME Team:

No! MEME uses a greedy algorithm. As soon as it has enough evidence to estimate the statistical significance of a motif candidate it will stop looking for further evidence for that motif, and start looking for other motifs. The statistical significance of a motif discovered by MEME is given by the E-value. An E-value of 1.3e-357 indicates that motif is highly statistically significant in your sequence data.

Keep in mind that MEME is performing de novo motif discovery without any reference to databases of known motifs, while CENTRIMO is looking for enrichment of motifs both from your MEME results and from databases of known motifs. You'll want to asses whether the highly significant motif reported by MEME is really some variant of the known motif reported by CENTRIMO. You can compare the logograms by eye, and then look at the TomTom results to see if they really represent the same motif.

> The known motif which was found with centrimo, indicates 1658 in "region matches" section. does it mean that 1658 of my sequences have the motif?

Clicking on the help button (the red question mark) for the "Region Matches" column of the CENTRIMO results provides the following text.

The number of (positive) sequences whose best match to the motif falls in the reported region.

Note: This number may be less than the number of (positive) sequences that have a best match in the region. The reason for this is that a sequence may have many matches that score equally best. If n matches have the best score in a sequence, 1/n is added to the appropriate bin for each match.

Furthermore, CENTRIMO is only considering motif matches that pass a score threshold. So, strictly speaking, CENTRIMO reporting 1658 regions matches means that 1658 of your sequences contain at least one match to the motif in their central region that passes the score threshold. Other sequences might also contain instances of the motif outside the central region, or with poorer matches to the motif position weight matrix than the score threshold will allow.

Look at the FIMO output included in the MEME-ChIP results to get a list of all motif matches in your sequence data.
Question: Hey, I have been using Gene Ontology to try to understand some differently express gene and I am now wondering if I should trust the GO annotation.

I was looking at serpinb2, for this gene I am getting GO term: GO:0005576 'extracellular region' and GO:0005615 'extracellular space'. However when I cross reference with Protein atlas, which have antibodies staining assays, then it seems to be clearly inside the cell around the golgi apparatus.

I have been trying to get the reference of those go annotation without success. Therefore which one should be trust, annotation or antibodies staining ?

Answer: I would encourage you to read about the Gene Ontology's evidence codes: http://www.geneontology.org/page/guide-go-evidence-codes

When you do your enrichment, you should see one of those codes beside your gene(s) of interest. The two GO terms that you've mentioned are very broad and undoubtedly have many 1000s of genes listed.

I just searched and SERPINB2 belongs to a few different GO categories; however, for both `GO:0005576 'extracellular region'` and `GO:0005615 'extracellular space'`, the evidence code is just 'TAS', i.e., '*Traceable Author Statement (TAS)*' (see through the link on evidence codes that I pasted above). Thus, you can make your own interpretation on that...

All of *SERPINB2*'s GO terms are here: http://amigo.geneontology.org/amigo/gene_product/UniProtKB:P05120

My general rule of thumb, and something that I tell all students that I train, is to never trust gene enrichment unless there is clear experimental evidence to back up the enrichment.

Hope that this helps.
Question: I have the following object in a basic DESeq2 bulk RNA-seq differential expression pipeline (human data). It filters out the genes that have low counts but on top of that I would like to remove a couple of genes that I know have issues in my dataset, and want to see how my analysis looks without them. I have the list of such genes in a vector named "genes" and it's encoded as gene symbols (I could transform them to EnsemblIDs if needed).

genesToRemove<- c("gene1","gene2","gene3","gene4","gene5","gene6")

dds <- DESeqDataSetFromHTSeqCount(sampleTable = mytable,directory = directory,design= ~ condition)

dds

class: DESeqDataSet
dim: 60725 326
metadata(1): version
assays(1): counts
rownames(60725): ENSG00000278625.1 ... ENSG00000277374.1
rowData names(0):
colnames(326): 9275 9351 ... 10146 10199
colData names(5): Condition Age ...

genes_to_keep <- rowSums(counts(dds)) >= 50
dds2 <- dds[genes_to_keep,]

I would like to do it at this point, after this code, so that then I keep going without them. The problem is that I am not sure how to access the part of the dds2 object where you have the genes in order to filter them out. Any thoughts? Thank you.

Answer: Hi!

In this case you could try this one:

#Obtain the indices of only desired genes
genesToRemove <- which(!rownames(dds) %in% genesToRemove)

#Cut your desired genes in the DESeq object
dds <- dds[genesToRemove, ]

#Verify that undesired genes are removed from DESeq object
genesToRemove %in% rownames(dds)

And the result must be `FALSE` for every undesired gene.

Best regards!
Question: Hi all,

I have a genomic regions sequenced from few accessions of a plant species. When I built a MSA of this region in different accessions. It looks like only the 4 genes present in this region are aligning. These four genes are less than 1 kb and the total genomic region is 16 kb. I know I have to provide more details to my query but does it look like something unusual (faster evolution) or it is normal that nongenic regions does not align in different accessions. If something is interesting then what I can do further.

Thansks,
Reetu

Answer: This is entirely normal from my experience with the A. thaliana genome of the Columbia cultivar. The following refers to A. thaliana BAC T10P11 I annotated in 1997.

Scroll down to gene T10P11.15, encoding LUMINIDEPENDENS. Look at the /note I added to the CDS for this gene. That reads, in part: 100% identical to peptide from cDNA U03456, from cultivar Wassilewskija (WS); exon 4 from cultivar Columbia is split into 2 exons in WS, while exons 11 and 12 in Columbia are joined as a single exon in WS.

The intron between exons 11 and 12 in Columbia is 135 bp in size, which is deleted in WS.

What to do further: Report these differences in sequence. See if the differences map very near markers for phenotypes indicative of the accessions you're studying (WS has darker green leaves than Columbia, if I recall correctly). Consider an analysis to incorporate these differences into a timeline of the emergence of the accessions. Lastly, look for evidence, like a footprint, of transposon activity - perhaps one was there and now left.
Question: I'm using TopGo to do GO analyses. When I get my enriched terms, the long term names get corrupted. ```I have some text then ...```. For example:

```
GO:0000122	negative regulation of transcription by ...
GO:0006886	intracellular protein transport
GO:0010976	positive regulation of neuron projection...
```

Any help would be appreciated?

My entire code:

```
library(topGO)
library(GO.db)
library(biomaRt)
library(Rgraphviz)
library(tidyverse)
library(dbplyr)

prebg_genes <- read.table("mart_mmu_protein_coding_genes.csv", header=F, sep=",")
bg_genes <- as.character(prebg_genes[,1])

pre_candidate_list_1 <- read.table("All_targets_ENSMUSG.txt", header=T, sep=",")
candidate_list_1= as.character(pre_candidate_list_1[,1])

db <- useMart('ENSEMBL_MART_ENSEMBL',dataset='mmusculus_gene_ensembl', host="www.ensembl.org")
go_ids <- getBM(attributes=c('go_id', 'ensembl_gene_id', 'namespace_1003'), filters='ensembl_gene_id', values=bg_genes, mart=db)
go_ids <- read.table("go_ids_mmu_protein_coding_genes.txt", header=T, sep="\t")
listAttributes(db)
View(go_ids)

gene_2_GO <- unstack(go_ids[,c(1,2)])

pre_keep1 <- candidate_list_1 %in% go_ids[,2]
keep1 <- which(pre_keep1==TRUE)
candidate_list_final1 <- candidate_list_1[keep1]

geneList1 <- factor(as.integer(bg_genes %in% candidate_list_final1))
names(geneList1) <- bg_genes

GOdata1 <- new('topGOdata', ontology='BP', allGenes = geneList1, annot = annFUN.gene2GO, gene2GO = gene_2_GO)

classic_fisher_result1 <- runTest(GOdata1, algorithm='classic', statistic='fisher')

weight_fisher_result1 <- runTest(GOdata1, algorithm='weight01', statistic='fisher')

applied to the topGOdata object.
allGO1 <- usedGO(GOdata1)
all_res1 <- GenTable(GOdata1, weightFisher=weight_fisher_result1, orderBy='weightFisher', topNodes=length(allGO1))

p.adj1 <- round(p.adjust(all_res1$weightFisher,method="BH"),digits = 4)

pre_all_res_final1=cbind(all_res1,p.adj1)
all_res_final1 <- pre_all_res_final1[order(pre_all_res_final1$p.adj1),]
write.table(all_res_final1,"summary_topGo_analysis_all_DE_miRs_targets.csv",sep=",",quote=F)
```

Answer: Crossposted on Bioconductor: https://support.bioconductor.org/p/131412/
Question: I have a large compressed VCF that needs to be compressed with bgzip so I can run bcftools. However, it isn't compressed with bgzip and if I decompress I will go over my quota on the cluster. Is it possible to compress a file without decompressing it fully and then compressing it again?

Answer: Assuming it's compressed with gzip, then zcat it to stdout which won't use up any of your disk quota and then bgzip.

```
zcat vcf.gz | bgzip -c > vcf.bgz
```