cannin/biostars_qa huggingface dataset:
Question: After peakcalling, i want to use igv to see my data(bedgraph). But i can not install IGV.
Do you have good suggestion?
Thanks!

Answer: If you are using Ubuntu 18.04 you can do `sudo apt install igv`

For previous versions you need something like:

sudo apt install default-jre
curl -LO http://data.broadinstitute.org/igv/projects/downloads/2.4/IGV_2.4.10.zip
unzip IGV_2.4.10.zip
cd IGV_2.4.10
sudo ln -s $PWD/igv.sh /usr/local/bin/
igv.sh&
Question: I have the code below to generate the heatmap with annotation. ![enter image description here][1]

library(pheatmap)
# Generate some data
test = matrix(rnorm(200), 20, 10)
test[1:10, seq(1, 10, 2)] = test[1:10, seq(1, 10, 2)] + 3
test[11:20, seq(2, 10, 2)] = test[11:20, seq(2, 10, 2)] + 2
test[15:20, seq(2, 10, 2)] = test[15:20, seq(2, 10, 2)] + 4
colnames(test) = paste("Test", 1:10, sep = "")
rownames(test) = paste("Gene", 1:20, sep = "")
# Add annotation as described above, and change the name of annotation
annotation <- data.frame(Var1 = factor(1:10 %% 2 == 0, labels = c("Exp1", "Exp2")))
rownames(annotation) <- colnames(test) # check out the row names of annotation
pheatmap(test, annotation = annotation)

I'm wondering whether there is a way to add the labels in the annotation bar in R code. Thanks in advance.
![enter image description here][2]

[1]: https://s22.postimg.cc/6beg2ma81/Picture1.png
[2]: https://s22.postimg.cc/z2b9su8dt/Picture2.png

Answer: use grid.text (part of base functions, earlier part of library "grid") (play with x and y coordinates):

pheatmap(test, annotation = annotation)
grid.text(levels(annotation$Var1), x=c(0.25,0.6),y=c(0.89,0.89), gp=gpar(fontsize=10))
Question: Hi everyone,

I was wondering if anyone is familiar of any annotation term in the human genome annotation e.g. from gencode or ensembl to be able to extract **intronless genes** and separate them from genes containing introns.

There is probably an automated way to extract intronless genes from the exon annotations in the gtf files. If there are any initial thoughts on this it would be much appreciated

Thanks in advance,

Sergio

Answer: A simple perl one-liner:

perl -lne 'if (/.+\texon\t.+gene_id "([^"]+)/) { $g{$1}++ }
END { foreach $i (sort keys %g) { print $i unless $g{$i} > 1 } }
' gencode.v26.annotation.gtf \
> intronless.txt

If you want to make sure that the code works, you can have both intermediate (exoncounts.txt) and final results (intronless.txt) to check manually the exon count:

perl -lne 'if (/.+\texon\t.+gene_id "([^"]+)/) { $g{$1}++ }
END { foreach $i (sort keys %g) { print "$i\t$g{$i}" } }
' gencode.v26.annotation.gtf | tee exoncounts.txt \
| perl -lane 'print $F[0] if $F[1]==1' > intronless.txt

Taking the same exon counting rationale into this `grep | cut | sort | uniq | perl` combination is even faster:

grep exon gencode.v26.annotation.gtf | cut -d'"' -f2 | sort | uniq -c \
| tee exoncounts.txt | perl -lane '$F[0] == 1 and print $F[1]' > intronless.txt

Important note: just realized that CDS lines do get through the previous grep, plus that HAVANA and ENSEMBL annotations may be redundant (therefore same exons could be counted twice), so the code should consider those issues in order to generate the proper output:

awk '{ if ($3=="exon") print $1, $4, $5, $10 }' gencode.v26.annotation.gtf | sort -u \
| cut -d'"' -f2 | sort | uniq -c | perl -lane '$F[0] == 1 and print $F[1]' > intronless.txt

Explanation: `awk` selects exon lines and prints only chromosome, start, end and geneid, `sort -u` collapses redundant exons, `cut -d'"' -f2` reduces output to geneids only, `sort | uniq -c` collapses same geneids while counting them, and `perl` prints geneids containing 1 exon only.
Question: I am trying to generate PON for WES samples following [GATK recommendations][1], they also have another explanation in this [Mutect2 article][2] but it's basically the same identical 3-steps procedure:

> ***step 1**. Run Mutect2 in tumor-only mode for each normal sample*:
`
gatk Mutect2 -R reference.fasta -I normal1.bam -max-mnp-distance 0
-O normal1.vcf.gz
`

> ***step 2**. Create a GenomicsDB from the normal Mutect2 calls*:
`
gatk GenomicsDBImport -R reference.fasta -L intervals.interval_list
--genomicsdb-workspace-path pon_db
-V normal1.vcf.gz
-V normal2.vcf.gz
-V normal3.vcf.gz
-V ...
`

> ***step 3**. Combine the normal calls using CreateSomaticPanelOfNormals*:
`
gatk CreateSomaticPanelOfNormals -R reference.fasta --germline-resource af-only-gnomad.vcf.gz
-V gendb://pon_db
-O pon.vcf.gz
`

I am using `gatk 4.1.7` (latest at the moment) but the output I got from step 2 (`GenomicsDBImport`) is a folder with some files in it, such as `vcfheader.vcf`, `vidmap.json` and what looks like a file for every chromosme with a `$` and contig boundaries specified in the BED file (e.g. `X$200786$155255277`).

If I try to pass this directory in the `-V` option of `CreateSomaticPanelOfNormals` (step 3 ) I got an error that the specified input is not a regular file, and GATK documentation confirms that `-V` is supposed to be a VCF file.

Does anybody, that maybe has generate PONs before or worked with this, knows what is the exact file output from step 2 that I am supposed to pass in step 3 `-V`?

Thank you very much in advance for any help!

[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360042479112-CreateSomaticPanelOfNormals-BETA-
[2]: https://gatk.broadinstitute.org/hc/en-us/articles/360035531132--How-to-Call-somatic-mutations-using-GATK4-Mutect2

Answer: I've done this quite recently, with what I hope is the latest version.
In my case I had to do this as below:

`-V gendb://pon_db`

the `gendb://` shouldn't be changed, and what follows is the path to the directory created by `GenomicsDBImport`
Question: How to Grep the complete sequences containing a specific motif in a fasta file using shell command? Also, I want to include the lines beginning with a `>` before these target sequences.

I found this post :https://www.biostars.org/p/274859/ similar to my problem but I'm looking for different motif:

My motifs look like that :

SXXXX(F/S)XXXL

All my fasta file in one line and I have more than 300 sequences:

for example:

my sequence :

>sp|Q9H257.2|CARD9_HUMAN RecName: Full=Caspase recruitment domain-containing protein 9; Short=hCARD9
MSDYENDDECWSVLEGSRVTLTSVIDRSRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9H37.2|CTYU_HUMAN
HHHSVLEGFRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9re7.2|CARer_HUMAN RecName
BKLSVLEGWRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

The result should be only the first two sequences because they have the motifs `SXXXX(F/S)XXXL`

>sp|Q9H257.2|CARD9_HUMAN RecName: Full=Caspase recruitment domain-containing protein 9; Short=hCARD9
MSDYENDDECWSVLEGSRVTLTSVIDRSRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

>sp|Q9H37.2|CTYU_HUMAN
HHHSVLEGFRVTLTSVIDRFRITPYLRQTKVLNPDDEEQVLSDPNLVIRKRKVGVLLDILQRTGHKGYVAFLESLELYYPQLYKKVTGKEPARVFSMIIDASGESGLTQLLMTEVMKLQKKVQDLTALLSSK

I tried these command but that returned all three sequences

grep 'S...F\|S\|L.\(.\)\1\{4\}' jara3.fasta -B 1 > jara4.fasta

Answer: I believe that {4} that you are using is an extended regular expression and you would need to either use egrep or the -E flag with grep.

I got it to work using this:

grep -E 'S[A-Z]{4}[FS][A-Z]{3}L' jara3.fasta > jara4.fasta

Hope that works!
Question: Hey
I already did automatic genome annotation using RAST software. I am looking for simple visualisation tool or anything which can provide me information about ORF total number of given dataset or would let me extract decription of all ORFs together.

It's pretty simple to do, but programs I am using are not counting orfs and they won't let me copy all ORFs description where I could just number it on word or any other tool. I am using genebank format.

I am not looking for tools like ORF FINDER, orfs are already predicted and translated in my dataset.

Answer: I solved this problem using program [SnapGene][1]. I put genebank sequences into the program. choosen chromosome -> features -> switch off "full descriptions" -> ctrl + a

You will have info about number of selected features. "(X - 1 - Y)/2" is number of your annotated ORFs.
X - number of features
Y - number of features that are not genes or CDS ( i. e regulatory regions) (to check if you have that type of features press "sort by gene" button

For smaller datasets you can just manually select all genes and program will count it for you, but it would be problematic with huge datasets.

I choose this pathway because I had my files in genebank format. For fasta files I would reccomend solution of Mr Dlakic, much quicker. Thank everyone for help

[1]: https://www.snapgene.com/
Question: Hi all,

I have a customed database and I used it to BLASTn against a bacterial genome. **I would like to extract the unmatched regions only.**

Is there a command line or another way to do it?

Thanks very much for your precious help!

Answer: 1. output in outfmt 6
2. convert to BED/GTF/GFF
3. `bedtools complement`
3. `bedtools getfasta`
Question: I did google this honestly, but I didn't understand what is the difference between all these names. I'd highly appreciate any help.

Answer: Conda is a package manager - a program where when you ask for a package to be installed, it will download an install the package you ask for and all of its depedencies (things it needs to run). The information on where to download things from and what the dependencies are for each package are stored in a database called a "channel".

bioconda is a conda channel, that contains the names, locations and dependencies of many bioinformatics tools. So when you say `conda install pysam -c bioconda` the conda program will go access the bioconda channel and ask 'where do I find the package "pysam", and what else needs to be installed?' and bioconda will reply "pysam is at https://blahblah.com/bioconda/pysam/versionXY, and when you download it, copy the files to direcoty ABC and then run script fgh.py. You will first need to install the packages htslib (version > 1.2.3) and samtools (version > 4.5.6) and python (version 3.8) and zlib (verion 9.10.11)". Conda then goes away and repeats that process for htslib, samtools, python and zlib until everything is ready to install, it downloads everything and installs it.

Anaconda is both a conda channel (that is a list of packages, where they can be downloaded from, and what their dependencies are) and an installable bundle that includes python, the most common data science packages for python (including things like scipy, numpy, pandas and matplotlib) and the conda pacakge manager.

miniconda is a different bundle that you can download from the same place as Anaconda, and it includes python and the conda package manager, but not the common data science packages, which you would have to install manually using the conda package manager if you later decided you wanted them.
Question: I have an alignment file for a protein encoded by about ~3000 bases (I have both bam and sam files), and I am looking to extract all the reads that from bases 877-981. Is there any way to do this with samtools?

I tried using this command:

samtools view 8.2alnCorrected.bam "877-981" > filteredreads.bam

However I got this message with an empty file:

region "877-981" specifies an unknown reference name. Continue anyway.

Answer: You need to specify the chromosome or contig as well. Your bam will have lines that look like:

chr1 12345 ...

So for reads are mapped to that chromosome or contig, you want:

samtools view 8.2alnCorrected.bam chr1:877-981
Question: Hello,

I have been using scanpy to analyze some single cell data using google colab.

Everything was working ok but all of a sudden I rerun the notebook and now when I try to import scanpy I get the error below. What could be causing this?

```
---------------------------------------------------------------------------
TypeError
Traceback (most recent call last)
in
2 import numpy as np
3 import pandas as pd
----> 4 import scanpy as sc
5 import matplotlib.pyplot as plt
6 get_ipython().run_line_magic('matplotlib', 'inline')

3 frames
/usr/local/lib/python3.8/dist-packages/scanpy/plotting/_utils.py in

33
34
---> 35 class _AxesSubplot(Axes, axes.SubplotBase, ABC):
36 """Intersection between Axes and SubplotBase: Has methods of both"""
37

TypeError: metaclass conflict: the metaclass of a derived class must be a (non-strict) subclass of the metaclasses of all its bases
```

Thank you

Answer: Hey y'all a fellow from scverse replied to my question and resolved our issue.

The issue seems to be that matplotlib version 3.7 is incompatible with scanpy version 1.9.1

As such, to resolve this issue you need to do is install a version of matplotlib thats older than 3.7. I used version 3.6 and that completely eradicated the metaclass issue.

pip install 'matplotlib == 3.6'

For further reference see the Github page:

https://github.com/scverse/scanpy/issues/2411
Question: After going through the manual I couldn't proceed with collapse replicates

sampleFiles <- list.files(path = "./", pattern = ".counts")
sampleNames <- gsub(".counts", "", sampleFiles)
sampleCondition <- c(rep("KO1", 4),rep("KO2", 4),rep("KO3", 4), rep("WT1", 4),rep("WT2", 4),rep("WT3", 4))
sampleTable <- data.frame(sampleName = sampleNames, fileName = sampleFiles, condition = sampleCondition)
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,
directory = directory,
design = ~ condition)
treatments <- c("KO1", "KO2", "KO3", "WT1", "WT2", "WT3")
library("DESeq2")
ddsHTSeq <- DESeqDataSetFromHTSeqCount(sampleTable = sampleTable,design = ~ condition)
colData(ddsHTSeq)$condition <- factor(colData(ddsHTSeq)$condition,
levels = treatments)
#Analysis using DESeq
dds <- DESeq(ddsHTSeq)
resultsNames(dds)
#Pre-filtering
dds <- dds[ rowSums(counts(dds)) > 1, ]
res <- results(dds)
#summarise some basic tallies
summary(res)

**This was the example from collapseReplicates, I am not sure where should I apply this step and how should I proceed with my dataset**

dds <- makeExampleDESeqDataSet(m=12)

# make data with two technical replicates for three samples
dds$sample <- factor(sample(paste0("sample",rep(1:9, c(2,1,1,2,1,1,2,1,1)))))
dds$run <- paste0("run",1:12)

ddsColl <- collapseReplicates(dds, dds$sample, dds$run)

# examine the colData and column names of the collapsed data
colData(ddsColl)
colnames(ddsColl)

# check that the sum of the counts for "sample1" is the same
# as the counts in the "sample1" column in ddsColl
matchFirstLevel <- dds$sample == levels(dds$sample)[1]
stopifnot(all(rowSums(counts(dds[,matchFirstLevel])) == counts(ddsColl[,1])))

Answer: In order to use `collapseReplicates`, you need a colData like this :

condition sample run
1 KO KO1 AdipoKO1a
2 KO KO1 AdipoKO1b
3 KO KO1 AdipoKO1c
4 KO KO1 AdipoKO12
5 KO KO2 AdipoKO2a
6 KO KO2 AdipoKO2b
...

Then you will collapse your "runs" (technical replicates) at the level of your samples (biological replicates) :

ddsColl <- collapseReplicates(ddsHTSeq, ddsHTSeq$sample, ddsHTSeq$run)
Question: Hi,
I need to create a table with how e-values are distributed for some sequences, as a way of reporting how conserved the sequences are.

I got some inconsistent results, and boiled it down to if I query a sequence alone, or if I query it together with other sequences. The result page has a drop-down menu where you can only pick a single query sequence. So I assume that it is independent of the other query sequences?

Here is an example to show it.
In the first test, I query ">1" alone, and the top hits are 4e-9.
In the second test, I query ">1" together with ">0" and ">2", and when I look at only ">2", the top hits are "3e-9"

First test set-up: ![https://i.imgur.com/E3xJqOj.png][1]

First test results: ![https://i.imgur.com/AOccFBj.png][2]

Second test set-up: ![https://imgur.com/uYYnnkC][3]
Second test results: ![https://imgur.com/v3SJUZ7][4]

I just did the same test with other sequences, and I got either 0.35 or 0.1 as the top hits. All settings are identical between the two searches. I just go to nucleotide BLAST, enter my queries, enter an organism, change to "blastn", and change the number of hits to 20.000. All other settings are the defaults.

So what is the **correct** way of doing this search? I'm so confused at the moment :<

[1]: https://i.imgur.com/E3xJqOj.png
[2]: https://i.imgur.com/AOccFBj.png
[3]: https://i.imgur.com/uYYnnkC.png
[4]: https://i.imgur.com/v3SJUZ7.png

Answer: Update. The NCBI help desk responded quickly, was able to replicate the bug, and quickly also addressed the bug.

For others that may have performed similar BLAST searches I paste their response below.

> The developers have addressed this issue.
>
> In summary:
>
>• The problem only occurred for BLASTN/megaBLAST searches.
>
>• It only happened if multiple queries were submitted at once. Results for the
> first query would be correct, but all other searches would use the
> search space for the first query instead of for each individual query.
>
> • It only affected the web page. Stand-alone BLAST+ does not have this
> issue.
>
>• Also are you aware of this paper:
> https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6662297/ This does have
> implications for E-values in some situations.
Question: So, first I tested what results I should get from the blastall program using the command line, with e-value 0.001:

C:\Niek\Test\blast-2.2.17\bin\blastall -p blastp -d C:\Niek\Test\arabidopsis-smallproteins.fasta -i C:\Niek\Test\arabidopsis-HD.fasta -e 0.001 -F F -m 8 -o C:\Niek\Test\arab-HD-smallproteins-notfiltered.out

and

C:\Niek\Test\blast-2.2.17\bin\blastall -p blastp -d C:\Niek\Test\arabidopsis-smallproteins.fasta -i C:\Niek\Test\arabidopsis-HD.fasta -e 0.001 -m 8 -o C:\Niek\Test\arab-HD-smallproteins-filtered.out

After that I made a local blast program, which works fine but it only found 91 results with e-value equal or lower than 0.001, where the results from the blastall via cmd gave around 140~ something results. I first thought it missed some, but all the e-values are different.

from Bio.Blast import NCBIStandalone
from Bio.Blast import NCBIXML

my_blast_db = r"C:\Niek\Test\arabidopsis-smallproteins.fasta"
my_blast_file = r"C:\Niek\Test\arabidopsis-HD.fasta"
my_blast_exe =r"C:\Niek\blast-2.2.17\bin\blastall.exe"

result_handle, error_handle = NCBIStandalone.blastall(my_blast_exe, "blastp",
my_blast_db, my_blast_file)

blast_records = NCBIXML.parse(result_handle)

E_VALUE_THRESH = 0.001
x = 0
for blast_record in blast_records:
blast_record = blast_records.next()
for alignment in blast_record.alignments:
for hsp in alignment.hsps:
if hsp.expect <= E_VALUE_THRESH:
print "==========Alignment========"
print "sequence:", alignment.title
print "length:", alignment.length
print "e value:", hsp.expect
x += 1

I first thought that the local blast from biopython uses a different algorithm, but at 'myblastexe =r"C:\Niek\blast-2.2.17\bin\blastall.exe"' I specify the same program but it should be the same. Then I thought it had something to do with the filtering option but I checked both filtered and unfiltered and wasn't any of that.

If you know why the local blast from biopython NCBIStandalone gives a different result than doing it directly in the cmd, please let me know.

Thanks in advance,
Niek

edit: I checked, and it seems that the NCBIStandalone filters out 100% identities, which the blastall called by cmd does not. However this doesn't explain why the e-values are so different.

Answer: One major problem is you are skipping half the results with this:

for blast_record in blast_records:
blast_record = blast_records.next()
...

It should be just:

for blast_record in blast_records:
....

Or if you would rather call the next method explicitly for some reason, something like:

while True:
blast_record = blast_records.next()
if blast_record is None: break
...

Secondary issue: blastall is now being phased out by the NCBI who call it "legacy" BLAST, they encourage people to use BLAST+ instead, in this case blastp at the command line. As a result, the old Biopython wrappers for calling "legacy" BLAST are all considered obsolete.
Question: Hello,

I have a table in R which consists of 70k rows and 37 columns. Lot of cells have **"./."** which I want to modify and make it as **"ab"** . I tried to use `gsub()` but it does not give me the required output.

I used :

file <- gsub("./.","ab",file)

I want the change to happen throughout the file.
Is there any other way with which I can modify it? Thanks in advance.

Input:
Eg:

S.no chr pos gene_name S1 S2 S3
1 1 1290 X ./. 1/1 ./.
2 1 5822 Y 0/1 ./. ./.

Output

S.no chr pos gene_name S1 S2 S3
1 1 1290 X ab 1/1 ab
2 1 5822 Y 0/1 ab ab

It can be either ab or NA

Answer: Try to use fixed match:

file <- gsub("./.","ab",file, fixed = TRUE)

Or

file[ file == "./." ] <- "ab"

-----------

**Edit:** Using example data provided by OP.

# example input data
df1 <- read.table(text = "
S.no chr pos gene_name S1 S2 S3
1 1 1290 X ./. 1/1 ./.
2 1 5822 Y 0/1 ./. ./.", header = TRUE, stringsAsFactors = FALSE)

df1
# S.no chr pos gene_name S1 S2 S3
# 1 1 1 1290 X ab 1/1 ab
# 2 2 1 5822 Y 0/1 ab ab

df1[, c("S1", "S2", "S3")][ df1[, c("S1", "S2", "S3")] == "./." ] <- "ab"
df1
# S.no chr pos gene_name S1 S2 S3
# 1 1 1 1290 X ab 1/1 ab
# 2 2 1 5822 Y 0/1 ab ab