cannin/biostars_qa huggingface dataset:
Question: Hi all!

I'm trying to do a heatmap with `plink`'s output for IBD calculation (`.genome`file).
This file has several columns, but only three are important to me, so my parsed input format looks like this:

IID1	IID2	PI_HAT
ID1	ID2	0.0163
ID1	ID3	0
ID1	ID4	0.0155
ID2	ID1	0.0096
ID2	ID3	0.0125
ID2 ID4 0.475
...

I would like to do a heatmap with the PIHAT values for all my IDs (I have hundreds). To my understanding, `R`asks for a matrix as input for the plot, but I'm not able to parse this input into the correct format (I actually get the heatmap, but all the values are wrong). Could someone please give me advice or if there's another way to do the plot, I'd be as well to try it. Thank you very much in advance!

Answer: Input data for `heatmap` should be numeric matrix, your data is long format data, we need to reshape it to wide, then convert to matrix for plotting, see example:

library(reshape2)

# example data
df1 <- read.table(text = "IID1 IID2 PI_HAT
ID1 ID2 0.0163
ID1 ID3 0
ID1 ID4 0.0155
ID2 ID1 0.0096
ID2 ID3 0.0125
ID2 ID4 0.475", header = TRUE, stringsAsFactors = FALSE)

#convert long-to-wide
x <- dcast(df1, IID1 ~ IID2, value.var = "PI_HAT")

# convert to matrix with column AND rownames
myM <- as.matrix(x[ , -1 ])
row.names(myM) <- x$IID1

# I am converting all NAs to 0, reconsider if this is suitable in your case.
myM[ is.na(myM) ] <- 0

#then plot
heatmap(myM)
Question: Hi

I have normalized a microarray data and matched the prob identifiers to gene symbol but now for some genes I have several probes. For example for gene A I have several matched probes. So I have repetition for gene A. How I can take mean over the expression of repeated genes please and having an unique value?

This is my expression matrix

> head(array[,1:10,1:5])
GSM482796 GSM482797 GSM482798 GSM482799
1 OR2T6 0.0171 -0.1100 -0.0394 -0.0141
2 EBF1 0.1890 0.0222 0.0832 0.0459
3 DKFZp686D0972 1.9400 0.2530 0.3770 0.8310
4 ATP8B4 -0.1490 0.0690 -0.0637 -0.0527
5 NOTCH2NL 0.1540 -0.3880 0.2160 -0.0812
6 SPIRE1 0.2920 0.1690 0.5500 0.1430

Answer: > limma::avereps(array,ID=rownames(array))
Question: Hi all,

I have bulk RNA-seq data with 12 samples - WT (x4), 'A' KO (x4), and 'B' KO (x4). I want to generate a 2D PCA plot (biplot) like below figure to look at the relationship between the samples.

![a 2D PCA plot][1]

I have tried an R package, 'PCAtools,' but it looks not work correctly as below.
![the PCA plot generated by the code below][2]

I have pasted my code and data below. I will very much appreciate it if you share any advice or suggestions.

Thanks in advance!

Joshua

library(PCAtools)
data <-read.csv("C:/.../all,gene_log2cpm_revised.csv", fileEncoding = 'UTF-8-BOM')

groups = c(rep("WT", 4), rep("A-KO", 4), rep("B-KO", 4))
cols = c('red', 'green', 'blue')[factor(groups)]

data$gene_name = as.numeric(as.factor(data$gene_name))

pca = prcomp(data)
pca$x
pca$sdev
biplot(pca, cex=0.7, scale=T, xlim=c(-0.6,+0.6))

Data file format:

![Data file format top genes][3]

[1]: /media/images/7893d708-b043-4fc5-a698-665c99df
[2]: /media/images/428ae629-d824-4fc0-a2a2-9c61f108
[3]: /media/images/a9d95265-adc6-444c-9aca-17b400a3

Answer: Hi Joshua,

Your second plot is not produced by my package, *PCAtools*. It is produced via base R functions.

Can you take a look through my vignette, please? - https://bioconductor.org/packages/release/bioc/vignettes/PCAtools/inst/doc/PCAtools.html#quick-start-deseq2

You probably need something like:

library(PCAtools)
data <- read.csv(
'C:/.../all,gene_log2cpm_revised.csv',
fileEncoding = 'UTF-8-BOM',
row.names = 1)

groups = c(rep("WT", 4), rep("A-KO", 4), rep("B-KO", 4))
cols = c('red', 'green', 'blue')[factor(groups)]

metadata <- data.frame(groups = groups)
rownames(metadata) <- colnames(data)
p <- pca(data, metadata = metadata)

biplot(p, colby = 'groups')

Kevin
Question: I want to recall SNPs from Illumina HumanOmni2.5-4v1. I have the raw data (`Grn.idat` and `Red.idat`) files, and also a matching `FinalReport.csv` which include the next columns:

SNP Name Sample Name GC Score Allele1- Forward Allele2- Forward Allele -Top Allele2-Top Allele1-Design Allele2-Design Allele1-AB Allele2-AB Theta R X Y X Raw Y Raw B Allele Freq, Log R Ratio

And I have thousands of such files. How can I get the number of calls for a specific SNP?

Should I use the raw idat files, or CSVs? And if the answer is the CSVs, then which column, and how to interpret it?

Thanks

Answer: I finally ended up using R package named crlmm. It has a function named genotype.Illumina that does the recalling given the raw idat files.
Question: Hi everyone. Many bioinformatics courses usually have the same items like matching aligments, genome assembling, sequence alignment, variant calling and so on. But, I am interested in algorithms on which for example BLAST or Artemis works (like choosing primer, finding splicing indexes, find secondary structures of a family of RNA molecules etc).

Could you provide me some books or other sources where I can read about these non-obvious algorithms.

Sorry if my question is silly, I know that I could google every algorithm I want to understand, but it would be better to have 2-3 sources with a good explanation.

Thanks a lot.

Answer: Basically, majority of the algorithms in bioinformatics can be found the pdf files as the following:

1. http://bioinformaticsalgorithms.com/
2. http://www.math-info.univ-paris5.fr/~lomn/Cours/BC/Publis/Complements/introductiontoBioinformaticsAlgorithms.pdf
3. https://mitpress.mit.edu/books/introduction-bioinformatics-algorithms
4. http://www.amazon.com/Introduction-Bioinformatics-Algorithms-Computational-Molecular/dp/0262101068
5. http://bix.ucsd.edu/bioalgorithms/
6. http://www.comp.nus.edu.sg/~ksung/algo_in_bioinfo/
Question: I am trying to analyze GSE 30321 datasets in R. However, it has 295 samples separated into two files.
I use

gset <- getGEO('GSE30321',GSEMatrix=TRUE)

In this way, gset has two elements. Does anyone know how to merge the two files into one and gset has only one element with 295 files? Thank you.

Answer: gselist = getGEO("GSE30321")
eset = combine(gselist[[1]],gselist[[2]])

This is usually all that is needed.
Question: Hello,

I tried to download the library ERX5671923 from SRA using fastq-dump with the --split-files option. It is a library from the Fly Cell Atlas (experiment ERP129698 in SRA).
I retrieved only one file per run (e.g ERR6032593).
As it is a paired end 10X V3 library I was expecting to retrieve 2 or 3 files (Read 1, Read 2 and potentially Index 1) but it contains only one single file with 91bp reads.
Do you have any idea if it is possible to use this file and, if yes, how to use it? I want to generate abundance matrix using kallisto/bustools.

Best Wishes,

Julien

Answer: That happens quite often that R1 is missing, don't ask me why. Good thing is that often submitters provide BAM files allowing reconstruction of fastq from there. That is the case [here][1] for all four accession numbers. See for example at the bottom of [here][2].

You can conveniently get bam files with `prefetch` from the sra-toolkit:

```bash
mamba install -c bioconda sra-tools
prefetch --type bam --max-size 9999999999 -O ./ ERR6032593
```

Sometimes in `Type` (see below) it doesn't say `bam` but something like `10X Genomics bam file`, for example [here][3]. Then you can use `--type TenX` with `prefetch` afaik.

![enter image description here][4]

Once you have the BAM files **and it is the BAM file from CellRanger** use the bam2fastq utility from 10x to convert the bam back to fastq: https://support.10xgenomics.com/docs/bamtofastq

If the BAM was made from alternative pipelines you will probably need to do custom parsing to recreate the R1 file as technically scRNA-seq (10x) is single-end sequencing using R2 while R1 is not used for the actual alignment but CB/UMI are processed differently. You probably need to access the tags that store CB and UMI sequences and recreate R1 accordingly, putting these sequences into the read positions where either CellRanger or your processing pipelines expect them. For example, 10x Chromium 3' v3 has CB in R1 position 1-16 and the UMI at 17-28, so that is relatively easy to parse from the BAM tags (I guess, untested, never done manually myself). But then again there are probably corner cases, so be careful.

See also: https://bioinformatics.stackexchange.com/a/15523

[1]: https://www.ncbi.nlm.nih.gov/Traces/study/?acc=ERX5671923&o=acc_s%3Aa
[2]: https://trace.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&acc=ERR6032593&display=data-access
[3]: https://trace.ncbi.nlm.nih.gov/Traces/index.html?view=run_browser&acc=SRR5167880&display=data-access
[4]: /media/images/6b04d15d-d005-4699-a986-3ce7c53f
Question: Hello !

I don't know what the difference between " ENSG00000002586.1" and "ENSG00000002586.19_PAR_Y" in Ensembl gene/transcripts ID means.

How do you deal with these ID with "_PAR_Y" in RNA-seq mapping process ?

Thanks,

Answer: Just found out that these suffixes are added in the GENCODE annotation files. We need to discuss how we deal with these because we can't be putting out IDs with one hand and not providing tools that can interpret them with the other hand.

I'm sorry for any confusion.
Question: Hi!

I am currently working to analyze alpha & beta diversity with microbiome data of 48 sample.
(Animal stool samples)

I finished working of Beta diversity and got the data of significance. (Weighted UniFrac, Bray-curtis)

However, I received this message and need some help from you.

The message is as below (From reviewer)
------
I highly recommend the use of **linear modeling (LM) or generalized linear modeling (GLM) which is** commonly used in microbiome studies **rather than a Wilcoxon rank-sum test**. This will allow you to better control for things that may impact your results such as the age animals. It is also important to control for site as this is known to significantly affect the gut microbiome and could be included as a random variable in a generalized linear mixed model (GLMM).
----------

I found that non-binomial regression or possion regression is the most commonly used glm for
microbial analysis and trying to use the models in comapring Alpha and Beta diversity difference of the 16s rRNA data.

However, my question is that what package in R or python is recommended when figuring out microbial
analysis. (Perhaps, many researchers in the metabarcoding field use some common statistical tools)

And, also I want to know that how can I modify (adding a dependent variable or etc. ) the data file
to use in the GLM model. (What variable should be added to compare in Alpha and Beta diversity analysis, and how can I calculate them?)

my metadata file is as below

![Meta data file to be used in GLM][1]

[1]: /media/images/478208ba-66c0-4776-a09d-8228c249

Any recommendation about using Package source, specific GLM model or the way to modify data -
would be a big help for me to solve the problem.

Thank you

Answer: In R, linear models are available with the function lm() and generalized linear models with glm(), both from package stats. For mixed effects, you can use the [lme4 package][1] (functions lmer() and glmer()). There's also the function glmnet() from the [glmnet package][2] if you want LASSO or elastic net regularization.

[1]: https://cran.r-project.org/package=lme4
[2]: https://cran.r-project.org/package=glmnet
Question: I have a list of genomic ranges mapped to hg19 . my data is in matrix format lets call it `ranges` which has 600,000 rows and 4 clumns
here is few row of my data

head(ranges)
chr start end strand
[1,] "chr1" "10025" "10525" "."
[2,] "chr1" "13252" "13752" "."
[3,] "chr1" "16019" "16519" "."
[4,] "chr1" "96376" "96876" "."
[5,] "chr1" "115440" "115940" "."
[6,] "chr1" "235393" "235893" "."
Is there a function that gets sequences and calculates GC content for each row( each range)
I would prefer that output be in a vector format
I would really appreciate your help

Answer: Hi,

You can use the `getSeq` and `alphabetFrequency` functions from Bioconductor's [Biostrings](https://bioconductor.org/packages/Biostrings/) package to obtain the sequence and compute the frequency of each nucleotide, and then, just divide the sum of G+C by the total number of nucleotides.

To use `getSeq` the easiest way is to create a GRanges object from your matrix. I use `toGRanges` from [regioneR](	https://bioconductor.org/packages/regioneR/) package, but there are other ways to build it.

library(regioneR)
library(BSgenome.Hsapiens.UCSC.hg19)

ranges <- matrix(c("chr1", "10025", "10525", ".", "chr1", "13252", "13752", "."),ncol = 4, byrow = TRUE)
colnames(ranges) <- c("chr", "start", "end", "strand")

#Build a GRanges from your matrix
ranges <- toGRanges(data.frame(ranges))

#Get the sequences and compute the GC content
freqs <- alphabetFrequency(getSeq(BSgenome.Hsapiens.UCSC.hg19, ranges))
gc <- (freqs[,'C'] + freqs[,'G'])/rowSums(freqs)

Here `gc` is a vector with the GC frequency of each range in the original matrix.

Hope this helps

Bernat
Question: The site that normally hosts the gmap/gsnap source code has gone away:

http://research-pub.gene.com

Does anyone know if this intentional or when it will be back?

Answer: This has now been resolved and http://research-pub.gene.com/gmap/src/ is available again.

For future reference raising an issue via this URL under the "Website Usability" seemed to get this to the correct people.

https://www.gene.com/contact-us/email-us
Question: Hi all,

Our lab has sequenced 20 different strains of the same bacterial species. I now want to analyse these data, i.e., link genomic data to virulence (some of the strains are virulent, others aren't) and other phenotype differences.

I started with de novo assemblies and mapping to an annotated reference genome. However, it is unclear to me how to proceed now in a good way. Most of the information for comparative genomics I find on the web applies only to the human genome or one-to-one comparisons of bacterial genomes.

Is there a standard workflow for this kind of analysis? Does anybody know of any good tutorials/references/articles? How would you proceed?

Answer: I don't think there is a single best standard approach to this, but [Bakker et al.][1] did a comparative genomics approach on Listeria with the goal of identifying virulence genes. You could simply try to replicate their methods on your data.

The main focus seems to be on the gene level, identifying a subset of virulence related candidate genes, and then look for whole gene deletions. Thereby I would also include biological knowledge on bacteria in general and knowledge about this specific species, as to what are good candidate sets of genes (e.g. Type III secretion system, etc.). The aim is to identify if there are differences in the inventory of genes.

In addition and in case, the whole gene deletion/insertion approach does not yield good candidate genes, I would run the whole raw data through a variant calling pipeline and look at smaller variation like small insertions or deletions, frame-shifts, mutations in CDS and promoter sequences. To accomplish this a standard pipeline (aka. [BWA][2] into [samtools][3]) could be used. You can then check variants common for certain phenotypes and check their effect on the protein level. Further, looking for non-coding elements might be interesting too.

[1]: http://www.biomedcentral.com/1471-2164/11/688/
[2]: http://bio-bwa.sourceforge.net/
[3]: http://samtools.sourceforge.net/
Question: Hi, I have data from a specific cell from mouse fed with a certain diet. I integrated 4 datasets that were measured at four different time for the integrated single-cell RNA seq analysis. I have been referring to the Seurat vignette : https://satijalab.org/seurat/v3.1/immune_alignment.html.

I am using SingleR to identify cell type for each cluster and I am wondering if I need to set `DefaultAssay` as "RNA" or "integrated". I tried both, but they gave me slightly different results for cell type identification.

Should I keep `DefaultAssay` as "RNA" or "integrated"?

Any thoughts and advice are greatly appreciated. Thank you.

Answer: The assay should be RNA, since SingleR expects expression values.
Question: Dear Friends,

I have a matrix like this:

Gene BRCA THYM TGHJ
ACC 23 21 7
XTG 12 13 9
CFG 45 4 8

The numbers are the snp count from vcf file.

I want to plot this in the form of a matrix with numbers colored based on their count; for example the highest number is colored in "Red" and then gradually the intensity of the color decreases with the decrease in number, so in this 45 is colored "Red" and 4 is colored with a very light color. Please let me know if am clear.

I am looking to plot this matrix using ggplot2, but other ways in R are also very welcome. A matrix like this:

![enter image description here][1]

[1]: https://s8.postimg.cc/65y00ijs5/Screen_Shot_2018-08-02_at_1.05.39_AM.png

Answer: I remember seeing a package for this exact task, but can't find it. So here is the starting steps, it will require more work to make it as pretty as you have in your example plot:

# example input
df1 <- read.table(text = "
Gene BRCA THYM TGHJ
ACC 23 21 7
XTG 12 13 9
CFG 45 4 8", header = TRUE)

library(ggplot2)
library(dplyr)
library(tidyr)

plotDat <- gather(df1, key = "Gene2", value = "value", -Gene)

ggplot(plotDat, aes(Gene, Gene2, col = value, fill = value, label = value)) +
geom_tile() +
geom_text(col = "black") +
theme_minimal() +
scale_fill_gradient2(low = "white", mid = "yellow", high = "red") +
scale_color_gradient2(low = "white", mid = "yellow", high = "red")

![enter image description here][1]

[1]: https://i.imgur.com/cIDojU0.jpg
Question: I'm interested in comparing the genes and genes sets regulated by two different drugs. I performed a microarray on Drug A vs control, and Drug B vs control. Using Limma, I calculated the differential expression for all 20,000 genes, and performed a pre-ranked GSEA. I now have ~700 gene sets that are depleted by Drug A, and ~800 gene sets that are depleted by Drug B (FDR<0.05). There are ~280 gene sets that are depleted by both Drug A and Drug B.

**Taking the pre-ranked GSEA output, how can I meaningfully compare these 280 gene sets depleted by both treatments to characterize similarities and differences between these treatments**.

If I am interested in PI3K signaling, could I compare the Normalized Enrichment Scores (NES) for all of the gene sets involved in PI3K signaling and perform a Wilcoxon test to determine if there is different enrichment between the two treatments?

Can I take the number of genes within a gene set regulated by each drug and meaningfully compare those gene sets?

Thank you for the help, a few pointers or a publication will set me on the right track. I've searched and can't find analyses like this.

Answer: I believe that gene set variation analysis (GSVA) will work here. This will allow you to compare two sets of enriched terms and pathways. It may take a day to learn the methods, but the tutorial is a good starting point.

Other than that, just a manual 'human' comparison by looking over the top enriched terms should suffice. At the end of the day, we have to learn to disengage the computer and start to interpret results ourselves.
Question: Hi all,
reading carefully the documentation of VarScan CNA pipeline I noticed in step 4 the following suggestion:

> If all of the data and segments are consistently above or below the
> neutral value (0.0), you can re-center the data points with VarScan
> copyCaller.

My data seem to belong in this category after plotting in R using DNAcopy package. All are consistently below 0.0.
So, my question is how should I know how much I do need to re-center my data?
I mean I can calculate this by eye, but how accurate can be that?
Is there any proper way to calculate that?

Thank you in advance.

Answer: Here's an old script I have laying around that does the recentering. Run it on the varscan copyCaller output, then repeat segmentation on the new recentered file.

https://gist.github.com/chrisamiller/f2c9e8bd565d500a8e8166fbd5eba2b1
Question: I'm trying to achieve what this post was looking for https://www.biostars.org/p/77012/

Currently this is my command:

bcftools mpileup -Ou --max-depth 8000 --min-MQ 30 --min-BQ 30 -f reference.fasta sample1.sorted.bam | bcftools call --ploidy 1 -Ou -mv | bcftools filter -s LowQual -e \'%QUAL<20\' > sample1.flt.vcf

#CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	sample1.sorted.bam
Imtechella_halotolerans_length_3113269	4051	.	C	T	41.4148	PASS	DP=2;VDB=0.02;SGB=-0.453602;MQ0F=0;AC=1;AN=1;DP4=0,0,0,2;MQ=42	GT:PL	1:71,0
Imtechella_halotolerans_length_3113269	4081	.	C	T	45.4146	PASS	DP=2;VDB=0.02;SGB=-0.453602;MQ0F=0;AC=1;AN=1;DP4=0,0,0,2;MQ=42	GT:PL	1:75,0

As far as I can understand the DP specified in the INFO is depth of coverage across all samples. How can I get the depth per sample info in the format/genotype field?

Any help appreciated!

Answer: Hi, you can specify this via the `mpileup` command. Please take a look at this snippet from my historical pipeline:

- ClinicalGradeDNAseq#L300-L316

That is:

bcftools mpileup \
--redo-BAQ \
--min-BQ 30 \
--per-sample-mF \
--annotate FORMAT/AD,FORMAT/ADF,FORMAT/ADR,FORMAT/DP,FORMAT/SP,INFO/AD,INFO/ADF,INFO/ADR \
-f "${Ref_FASTA}" \
Aligned_Sorted_PCRDuped_FiltMAPQ.bam |\
bcftools call \
--multiallelic-caller \
--variants-only \
-Ob > Aligned_Sorted_PCRDuped_FiltMAPQ.bcf ;

Kevin
Question: Is there a software package for GWASs that can include an interaction effect (SNP x predictor) in the regression using dosage data?

Answer: Can you use [PLINK2](https://www.cog-genomics.org/plink2/assoc#linear)?

That provides an interaction term for the linear regression... In terms of the dosage data, do you have to use the probabilistic data, or would best-guess hard-calling work?
Question: Hi everyone, I was wondering if someone can point me to a good explanation where i can better understand how the reads in my fastq files are created.

in the fastqc results from the data set I got, we assume that a lot of the adapter was sequenced ([2]. image). I'm not talking about read-thourgh, but the adapter somehow was sequenced. At least we think this is what we have.

I have tried to understand in which direction the reads in my fastq files are read. ![example Overrep. sequences ([1]. image). What I don;t get is on which side do I find my adapters, if any are left?

Are they at the beginning of my read so should I crop head of the read, or should I trim the end of the of the read, as the adapters are there?

Is there a good explanation for that somewhere?

thanks

george

![1]
![2]

[1]: /media/images/dd7a724b-980d-4bd9-8291-5ee57fc3
[2]: https://i.postimg.cc/9MZhcLMj/per-base-sequence-content.png

Answer: In Illumina sequencing adapter are always going to be present at 3'-end of the read (unless you are using some modification of standard procedure). You can also have adapter dimers (i.e no insert). That said, what kind of data is this. Looks like you can almost read the sequence looking at the plot so amplicons perhaps?