cannin/biostars_qa huggingface dataset:
Question: Hello all,

I've posted the question in Stackoverflow but I thought I might get more responses here.

I was able to load my csv file into a numpy array:

data = np.genfromtxt('csv_file', dtype=None, delimiter=',')

Now I would like to generate a heatmap.

I have 19 categories from 11 samples, along these lines:

COG station1 station2 station3 station4
COG0001 0.019393497 0.183122497 0.089911227 0.283250444 0.074110521
COG0002 0.044632051 0.019118032 0.034625785 0.069892277 0.034073709
COG0003 0.033066112 0 0 0 0
COG0004 0.115086472 0.098805295 0.148167492 0.040019101 0.043982814
COG0005 0.064613057 0.03924007 0.105262559 0.076839235 0.031070155
COG0006 0.079920475 0.188586049 0.123607421 0.27101229 0.274806929
COG0007 0.051727492 0.066311584 0.080655401 0.027024185 0.059156417
COG0008 0.126254841 0.108478559 0.139106704 0.056430812 0.099823028

I wanted to use matplotlib colormesh.

all the examples I could find used random number arrays.

I can get the plot easily with random numbers, however I can't get my csv file to plot. first it refuses to reshape. I have NaNs there so I tried masking but that failed too. Also, I had to delete the header and first column, is there a way to leave them and get labels for the axes? I've edited the original question to include an excerpt of the csv file.

any help and insights would be greatly appreciated.

many thanks

Answer: Here's a nickel, kid, go get yourself a better plotting library

> library(ggplot2)
> foo = read.table('foo.txt', header=T)
> foomelt = melt(foo)
Using COG as id variables
> ggplot(foomelt, aes(x=COG, y=variable, fill=value)) + geom_tile() + scale_fill_gradient(low='white', high='steelblue')
> ggsave('biostar.png')
Saving 7.97" x 7.75" image

ggplot2 is plotting heaven and *way* better than matplotlib. Use rpy2 to run from python - they even have ggplot2 examples in the docs.

![](https://imgur.com/fGTZT.png)
Question: I have made a package to which I call the QuasR package among others, therefore, I have put it both in the DESCRIPTION file and in the NAMESPACE file. The DESCRIPTION file is as follows:

> **LazyData**: true
>
>**Imports**: Rcpp,
Biostrings,
bedr,
devtools,
tidyverse,
dplyr,
phangorn,
ggseqlogo,
metan,
ggpubr,
scales,
ggplot2
>
> **Enhances**: parallel
>
> **Remotes**: bioc::release/Rbowtie,bioc::release/QuasR,bioc::release/GenomicAlignments,bioc::release/Rhtslib, bioc::release/rtracklayer,bioc::release/ShortRead, bioc::release/BSgenome,bioc::release/GenomicFeatures, bioc::release/VariantAnnotation,bioc::release/GenomicFiles,bioc::release/dada2, bioc::release/QuasR
>
> **LinkingTo**: Rcpp
>
> **RoxygenNote**: 7.2.0

The problem is that when I Dokerize this package, the following error occurs:

`Error in loadNamespace(i, c(lib.loc, .libPaths()), versionCheck = vI[[i]]):
there is no package called ‘QuasR’`

How can I call this package in the DESCRIPTION file?
I tried to insert as Imported package, with the same result...

Answer: Try to add a line with `biocViews:` before `Imports:` as suggested in this [answer][1]. Then you simply add the bioconductor libraries as `Imports:`

[1]: https://bioinformatics.stackexchange.com/questions/3365/r-package-development-how-does-one-automatically-install-bioconductor-packages
Question: Hi everyone,

I have used the "findOverlaps" function in R to find which positions of my two datasets overlap. I have also used "countOverlaps" to see how many overlaps I have.

What I want to do now, is to find the exact coordinates which overlapped. I was looking how to do it using findOverlaps but it has no such options and by doing a google search, I couldn't find much help.

A little help on this would be greatly appreciated.

Thank you

Answer: I have done this with dummy data.

Lets say we have we have regions A (groupA) and regions B (groupB).

groupA <- data.table(
chr = rep("chr1", 5),
start = seq(0, 1000, 100),
end = seq(50, 1050, 100))
setkey(groupA, chr, start, end)

groupB <- data.table(
chr = rep("chr1", 5),
start = seq(25, 1025, 100),
end = seq(75, 1075, 100))
setkey(groupB, chr, start, end)

Check:

1. If your datasets are data.table `class(groupA)`, if not do `setDT(groupA)`
2. If keys are chr start end, if not do `setkey(groupA, chr, start, end)`

# Find overlaps
over <- foverlaps(groupA, groupB, nomatch = 0)
# Extract exact regions
over2 <- data.table(
chr = over$chr,
start = over[, ifelse(start > i.start, start, i.start)],
end = over[, ifelse(end < i.end, end, i.end)])
Question: I am running DESeq2 to find DEGs between multiple samples, but I'm not able to decide what type of design to use, and how to arrange my data?

My data the following categories-

----------
**1.DISEASE SUBTYPE | 2. TYPE OF MUTATION**

A | mut1 | mut2 | mut3

B | mut1 | mut2 | mut3

C | mut1 | mut2 | mut3

----------
There are three different mutation backgrounds (mut1, mut2, mut3) common in each of the disease subtype (A, B, C).

**I want to compare my NORMAL HEALTHY samples with each of the Amut1, Amut2, Amut3, Bmut1, Bmut2... and so on. (and also inter-category comparisons)**

How should I arrange/tidy my data for DESeq2, and what should I write for the design?
Should I just compare them pairwise separately, or make a complex design for R?

Any kind of help is appreciated, I'm just a beginner. Thank you!!!

Answer: I would try out two approaches.

The first approach would be combining your disease subtypes and mutation factor levels for each sample. Your sample sheet would look as follows.

condition
sample_1 A_mut1
sample_2 A_mut2
sample_3 A_mut3
... ...

The regression formula would then just be `~ condition`.

The next approach uses interaction terms in the regression, so would model both the main effects of each mutation and each disease subtype, as well as the differences in effect for each disease subtype based on the mutation. Your sample sheet would now look like the following.

subtype mutation
sample_1 A mut1
sample_2 A mut2
sample_3 A mut3
... ... ...

For this analysis the regression formula would be `~ subtype + mutation + subtype:mutation`.

More information on multi-factor design can be found in the [DESeq2 documentation](https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html#multi-factor-designs), and a detailed explanation for working with interaction terms can be found in the help documentation for the results function, `help("results", package="DESeq2")`.
Question: I am new to analysing ATAC-seq data.

As mentioned on https://www.biostars.org/p/209592/, there seem to be two ways to use MACS to analyse ATAC-seq data.

1. Utilising the --shift -100 --extsize 200 command
This, I believe, is to find where the cutting sites are.

2. Utilising the --shift 37 --extsize 73 command
This is to find the nucleosomes since the DNA is wrapped around nucleosomes is circa 147bp.

I specifically want to document regions of open chromatin, particularly enhancer regions.

I have read in the literature that open chromatin regions will give reads <100bp. If this is correct, do I have to filter my bam file to look at only this length.

Surely if I map all of the my data, not filtering for read size I will end up with genome-wide ATAC-coverage? Looking at my bam file following mapping to the genome, my read sizes range from 50bp-450bp (with a peak at around 70bp).

Answer: All these macs2 options that people use, combining shift and extsizes in my experience do not make a notable difference. As spacemorrissey mentioned, unlike ChIP-seq (with typical short single-end reads) you want to scan ATAC-seq data for cutting sites within the accessable chromatin. I always use `macs2 callpeak -g hs --nomodel -f BAMPE`, which skips the macs-typical shifting model and piles up the entire paired-end fragment size.
For reliable enhancer calling, ATAC-seq is the wrong assay anyway, because even though active enhancers are accessible, by far not every accessible region is an enhancer (promoters, silencers, insulators etc...). What exactly is the question you are working on?
Question: Hello,

I am trying to recapitulate a figure in the following paper; figure 4F (both left and right panels):

https://www.cell.com/cell-reports/pdf/S2211-1247(21)01473-X.pdf

Does anyone know how to do these types of diagrams? I'm assuming they're done in Cytoscape. My data includes TFs in 4 different cell types, as well as the target genes (TFs) of those TFs in all 4 cell types, and then a column designating if that interaction if activating or repressing.

Currently I can create a circle network diagram in Cytoscape, but the TFs aren't organized by cell type other than the node color, so it just looks like one big circle. It also filtered any repeat TFs to just one node, so if I have the same TF in multiple cell types I can't visualize that...

Thank you!

Answer: Yes. This figure was made using Cytoscape (it is mentioned in the Methods). Here is how I would do it.

Data prep:
As you noted, you'll need unique node names for every TF you want represent. This means concatenating gene names with cell types in your case, e.g., Sox9_RPC and Sox9_MG. So, your data would look something like this:

Sox9_RPC | activates | Nfia_MG
Sox9_RPC | represses | Pax6_TT
Nfib_MG | activates | Lhx2_TT

In Cytoscape:

1. You can load these data into Cytoscape using File>Import>Network from file.
2. Next, you want to select the nodes for each of the cell types. Using the search bar, you can use `name:*_RPC` to select all the nodes with "_RPC", for example.
3. Then run the Circular Layout *on selected nodes only"*. Move that circle off to the side while it is selected. Repeat the selection and layout of each of the cell types.
4. You can use the Style panel to map edge color to the interaction type (a discrete mapping for "activates" or "represses"). You can also use bypass node border color to set the color of nodes per cell type when selected.
5. You can use the Annotation panel to add large font labels like "RPC" and "MG" for your cell types on the network (like in the figure).

These steps should reproduce that figure.
Question: Hi Biostars,

I have been using trimmomatic for quite some time, but realized that I don't get something.
So one can use LEADING and TRAILING options to remove bases from the beginning and end of the read, respectively.
From manual:
LEADING Remove low quality bases from the beginning. As long as a base has a value below this
threshold the base is removed and the next base will be investigated.

My question is:
"From the beginning" means from the beginning until the end of the read? If yes, then what is the meaning of having TRAILING option if the whole read is scanned?
Otherwise, until which base does trimmomatic scan by LEADING option?

Cheers,

Answer: According to the [documentation][1] both both options take the argument "quality":

LEADING:quality

leading: Cut bases off the start of a read, if below a threshold quality

quality: Specifies the minimum quality required to keep a base. Remove low quality or N bases.

Also from [the manual][2]: LEADING - Remove low quality bases from the beginning. As long as a base has a value below this threshold the base is removed and the next base will be investigated.

"LEADING 3" would delete all bases below a quality threshold of 3 or that are N, beginning at the first base and continuing until the first base that is at least a quality of 3 and is not N

AAAGGGTTT
012345678 - Leading 3 would cause the deletion of AAA

AAANNNTTT
012345678 - Leading 3 would cause the deletion of AAANNN

AAAGGGTTT
123456789 - Leading 3 would cause the deletion of AA

Same happens with trailing but from the other end

[1]: http://www.usadellab.org/cms/?page=trimmomatic
[2]: http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/TrimmomaticManual_V0.32.pdf
Question: I am working on some human WGS data hg38 assembly. I was wondering: is there any BED format to exclude low-quality sites?

Thank you very much in advance!

Answer: the complement of wgs_calling_regions.hg38.interval_list under https://console.cloud.google.com/storage/browser/genomics-public-data/resources/broad/hg38/v0/ ?
Question: I have a FASTQ file and I'm able to run the FASTQC program to analyse the file. but when I use trim_galore, FASTQC (or the FASTQC option in trim_galore) is not working anymore.

$ fastqc ./sub1_val_1.fq.gz

This is the output:

Started analysis of sub1_val_1.fq.gz
Analysis complete for sub1_val_1.fq.gz
Failed to process file sub1_val_1.fq.gz
java.lang.ArrayIndexOutOfBoundsException: -1
at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.calculateDistribution(SequenceLengthDistribution.java:100)
at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.raisesError(SequenceLengthDistribution.java:184)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.startDocument(HTMLReportArchive.java:336)
at uk.ac.babraham.FastQC.Report.HTMLReportArchive.(HTMLReportArchive.java:84)
at uk.ac.babraham.FastQC.Analysis.OfflineRunner.analysisComplete(OfflineRunner.java:155)
at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:110)
at java.lang.Thread.run(Thread.java:695)

Is the Failed to process file an error because the version is not correct between trim_galore and FastQC?

I found [this][1], [but that wasn't that helpful][2].

I'm using FastQC v0.11.5 and trim_galore v0.4.1.

I subsetted a library (reads in paired-end) using this:

seqtk sample -s100 ./SRR2937435_1.fastq.gz 10000 | gzip > sub1.fastq.gz
seqtk sample -s100 ./SRR2937435_2.fastq.gz 10000 | gzip > sub2.fastq.gz

The sub1_val_1.fq.gz file was after passing sub1.fastq.gz into trim_galore. FastQC with sub1.fastq.gz is working.

[1]: http://seqanswers.com/forums/archive/index.php/t-4846.html
[2]: https://bugs.launchpad.net/ubuntu/+source/fastqc/+bug/1443275

Answer: I found the answer: You have to uncompress it. Probably, trim_galore is only working with tar.gz and not fastq.gz.

gzip -d -k sub1.fastq.gz > sub1.fastq
y # to accept to overwrite
gzip -d -k sub2.fastq.gz > sub2.fastq
y # to accept to overwrite

trim_galore --illumina --paired --fastqc sub1.fastq sub2.fastq
Question: Hello, I'm fighting with my awk command since yesterday.

I have a file (locus.txt), this is some IgH locus from mm10 (I don't have header but you have : chr, start, end, strand and name_of_the_locus, separated by tab)

chr12	113363298	113365156	-	gamma3
chr12	113330756	113338695	-	gamma1
chr12	113308036	113314227	-	gamma2b
chr12	113274557	113277035	-	gammaepsilon
chr12	113260153	113264625	-	alpha
chr12	113289248	113295541	-	gamma2a
chr12	113423027	113426701	-	muIgh
chr12	113225832	113255223	-	3'RR
chr12	113416247	113418358	-	IgD

What I want to do is to grab the minimum position in this file, so the minimum position in start column (second column : `113225832`, for 3'RR)

Then, I want to substract all my position with this minimum and rearrange the file like this

gamma3	137466 139324
gamma1	104924 112863
...etc

----------

**What I have tried so far**

Search for minimum value, saved in $min :

min=`awk -v min=1000000000 '{if($2

Answer: To pass a variable to awk you can use the `-v` option like this (not tested):

awk -v min=$min -F $'\t' '{$1=$4=""; print $5 "\t" $2 - min "\t" $3 - min}' locus.txt
Question: On the Julia forum, there was the suggestion of creating a new "flavor" of the BioStar Handbook, but using Julia as programming language for the code examples.

https://discourse.julialang.org/t/biostar-handbook-computational-genomics-and-julia-to-be-or-not-to-be/25732

Is this allowed under the current license, and if so, would there be any interest from the current authors in collaborating on this?

Answer: Code examples may be translated in any other language,

What is not permitted is re-publishing another version of the book where any of the content is copied verbatim and only the code is changed.

As always *Fair Use* applies and I consider myself a supporter of [Fair Use](https://en.wikipedia.org/wiki/Fair_use) policies.
Question: I'm going through the instructions page on https://gatkforums.broadinstitute.org/gatk/discussion/1601/how-can-i-prepare-a-fasta-file-to-use-as-reference

Specifically, the command I don't see how to do is:

java -jar CreateSequenceDictionary.jar R= Homo_sapiens_assembly18.fasta O= Homo_sapiens_assembly18.dict
[Fri Jun 19 14:09:11 EDT 2009] net.sf.picard.sam.CreateSequenceDictionary R= Homo_sapiens_assembly18.fasta O= Homo_sapiens_assembly18.dict
[Fri Jun 19 14:09:58 EDT 2009] net.sf.picard.sam.CreateSequenceDictionary done.
Runtime.totalMemory()=2112487424
44.922u 2.308s 0:47.09 100.2% 0+0k 0+0io 2pf+0w

I think that `CreateSequenceDictionary.jar` comes from Picard, so I downloaded that from https://broadinstitute.github.io/picard/, but I don't see `CreateSequenceDictionary.jar` anywhere in the directory. However, I do see

`CreateSequenceDictionary.java`

I assume that `.jar` files are analogous to C executables, and `.java` files are analogous to `.c` human-readable code.

Going through the Picard readme file, I see that I should execute `./gradlew shadowJar` but this build fails on two different computers that I'm on.

So I can't make/get `CreateSequenceDictionary.jar`

I'm at a loss, how do I generate this dict file?

Answer: Picard is a wrapper command that will run the subcommands. There is no single executable for each subcommand. Simply run `picard.jar` and then check the output printed to screen for the subcommand you want. Then run `java -jar picard.jar `
Question: where to get the file of hg19 exon, intron, UTR region ?

I read lots of post in biostar. However it's been a long time post.

Answer: Get the respective GTF (annotation) file for your genome. Once you have this you can basically follow https://www.biostars.org/p/112251/#314840 to get the respective features. GTFs can be found at NCBI, Ensembl or GENCODE.
Question: Hi,

my question is regarding the most common diff. expr. tools used (edger, deseq2, limma). As they are created before UMIs existed (where the overall set of data includes quite low copy numbers on all entries which seem like rare low expressed genes) did such data can be used with their statistical methods and normalization and get proper DEG results? I'm not a statistician but know that changing the scale of the data, dispersion, etc may cause bias in the results?

Answer: UMI counts will work perfectly fine for these tools and should actually better match the statistical assumptions underlying them than regular read counts (due to the additional noise from PCR duplicates here). The issue is more when you get to drop-outs (missed signal), which can occur in things like scRNA-seq. In practice that still seems to be handled well enough by these tools.
Question: Hi Biostars,

Is it fine to use STAR for bacterial data (no splicing)? Any comments/suggestions are highly appreciated.

Thanks

Answer: It should work but you need to specify `--alignIntronMax 1` to force STAR to avoid splice alignments.

Also during the genome index generation step you should put `--genomeSAindexNbases` to `min(14, log2(GenomeLength)/2 - 1)` if your bacterial genome of interest is relatively small. GenomeLength is in base pairs
Question: Dear all,

considering a RNA-seq experiment and analysis that provides the expression values as TPM, please would you let me know what is a minimum TPM value in order to consider a gene to be expressed ?

talking about RPKM.FPKM units, I remember that a gene was considered expressed if RPKM (or FPKM) > 1 ... thanks a lot,

-- bogdan

Answer: I do not believe there is any definitive answer. There are so many factors that go into each experiment such that it is difficult to pick a value. A RPKM / FPKM value of 1 seems quite low to me, i.e., in 'error' territory. What you have to consider is the distribution of your data and the suitability of it for whatever downstream tools you will use. If including low-count / low expressed genes is going to distort your data distribution and introduce biases, then you need to remove them - check via histograms.

From RNA-seq, most genes are lowly expressed, possibly due to transcriptional 'noise' more than anything else. I say 'noise' in the knowing that they may reflect genuine transcription but have no regulatory function and are artifacts of other transcriptional processes that have occurred. They may also reflect regions where TF binding and/or promoter activity was weak.

So, you have the liberty to choose your own cut-off for TPM and state it in the methods. :)

Please take the time to read Gordon's answer, here: https://support.bioconductor.org/p/98820/#98875

Kevin

Edit: another interesting discussion: https://www.researchgate.net/post/How_to_determine_whether_a_gene_is_expressed_in_RNA-seq2
Question: There exists a lot of literature on distinguishing driver mutations from passengers. I am trying to build my own deep learning model to do the same. I am facing some potential issues. First, I have downloaded COSMIC mutation data and used the FATHMM labels to designate drivers (positive examples) in my dataset. I am sceptical to use passengers from COSMIC as they may be false negatives. So I turned to the 1000 genome project to download SNVs (to construct my negative examples). I am unsure if this is correct, however, I have seen some papers do the same. Do I need to apply any filters on the 1000 genome SNV data to construct the final dataset? One such paper talks of using SNVs with a global minor allele frequency≤1%.

Answer: This area of research is difficult because each person has their own opinion about what constitutes a 'driver' and 'passenger' mutation. For example, I disagree with the 1000 Genomes approach because I already know that some 1000 Genomes polymorphisms that have appreciable minor allele frequencies of around 15% in Caucasians can drive ER-positive breast cancer. To be frank: we don't have a clue what >90% of the variants listed in dbSNP / 1000 Genomes are doing. A large proportion of them could ultimately be driving cancer and other diseases.

If I were you, I would not try to define on my own what is / is not a driver and passenger. Why not utilise the work that has already been published and then follow their guidance for your deep learning model?

Look at this paper, published in the highly reputable *Cell* journal:

- Comprehensive
Characterization of Cancer Driver Genes and Mutations.

Here is another one, published a couple of weeks ago in *Nature Genetics*:

- Identification
of cancer driver genes based on nucleotide context

I also just stumbled upon this online platform, which focuses on drivers:

- intOGen

Build upon the work that is already out there. Then, at least, if you try to publish your work, it will be more difficult for reviewers to criticise you.

By the way, I am not sure why you are using FATHMM, or did you mean FATHMM-MKL, as mentioned HERE on COSMIC's page. GWAVA and Funseq2 were designed for somatic mutations ( see here: https://www.biostars.org/p/286364/#286483 ).

Kevin