cannin/biostars_qa huggingface dataset:
Question: I just start learning bioinformatics from RNA-seq/DEGs. I am trying to replicate an article. I very much appreciate any help. I used salmon in combination with tximport and DESeq2. PCA plot showed 77%PC1 8%PC2. Apparently, different conditions are mixed and not well separated. Is it healthy to continue this analysis? ![PCA][1]

library(tidyverse)
library(tximport)
library(DESeq2)

sample_table <- read.csv("SraTable.txt") %>% arrange(sample_table$Run) %>% select(Run)

sample_table$Run <- c("pre_ischemia1","ischemia1","pre_ischemia2","ischemia2","pre_ischemia3", "ischemia3","pre_ischemia4","ischemia4","pre_ischemia5","ischemia5")

SRR10700832_quant/quant.sf sample_files <- paste0(c("SRR10700832_quant","SRR10700833_quant","SRR10700835_quant","SRR10700836_quant","SRR10700838_quant", "SRR10700839_quant","SRR10700841_quant","SRR10700842_quant","SRR10700844_quant","SRR10700846_quant"), "/quant.sf")

names(sample_files) <- pull(sample_table, Run)
gene_map <- read.csv("gene_map.csv", col.names = c("enstid", "ensgid"))
count_data <- tximport(files = sample_files, type="salmon", tx2gene= gene_map, ignoreTxVersion = TRUE)
sample_table <- as.data.frame(sample_table)
sample_table$condition <- factor(rep(c("pre_ischemia", "ischemia"), times=5), levels = c("pre_ischemia", "ischemia"))

deseq_dataset <- DESeqDataSetFromTximport(txi= count_data, colData = sample_table, design = ~condition)
deseq_dataset <- DESeq(deseq_dataset)
vst <- varianceStabilizingTransformation(deseq_dataset)

plotPCA(vst, intgroup='condition')
results <- results(deseq_dataset)

summary(results)

`res <- results(deseq_dataset)` yielded following;

out of 38464 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up) : 1, 0.0026%
LFC < 0 (down) : 0, 0%
outliers [1] : 263, 0.68%
low counts [2] : 0, 0%
(mean count < 0)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

There is only 1 gene upregulated, zero down and zero low counts which are not likely. Which step could be the problem ? I would like to know sanity checkpoint in the entire flow.

Thank you for your time

P.S sorry it this hurts your eyes. I couldnt format the code in the post.

[1]: /media/images/2c5350ce-6c6e-4c8e-afb7-8034d5b2

Answer: Hello,

I just played with these data a bit as the salmon quantification directories were available at GEO, so this was no effort.

It seems the authors filtered on p-value based on the methods text...

The significantly differentially expressed genes (DEGs) between groups were defined at cut-off criteria of |log2 fold-change| ≥ 1 and p-value < 0.05.

...instead of padj. When you check some of the top20 DEGs that they provide in Supp. Table S3 you will see that these "significant" ones only pass the p-value cutoff but all have a padj >>> 0.05, so there indeed simply are no DEGs at standard cutoffs, at least not without correcting for potential batch effects. I did not look into the paper at any detail, but this at least seems to "replicate" the results they have. Not discussing now whether this strategy they chose was good or not.
Question: Hi,

I would like to know whether one needs to sort and to index the .bam files again after removing the duplicates in Samtools using `samtools rmdup -S`.

I am under the impression that since we have already sorted and indexed the .bam files, and we are just removing the duplicates, there would be no need to do the sorting and indexing again. Am I wrong?

Thanks in advance for your thoughts!

Answer: "sort"

no: your bam is already sorted.

"index"

yes if you need to access some specific genomic regions in your bam. The previous index of the bam cannot be used with the new rmdup.bam
Question: Hi!

I have assembled my PacBio FASTQ reads with Canu. Now I would like to polish\correct the assembly by mapping these PacBio FASTQ reads on the assembly itself. I heard about Quiver in the SMRTanalysis pack, but I'm wondering if there is any alternative software.

Thanking you in advance for your help!

Answer: Canu uses alignments between pacbio reads to correct away most sequencing errors in raw reads. You can't do much better without using the signal information. To achieve high-quality consensus, Quiver is your only choice.
Question: Does anyone have a good idea to combine a data frame with one containing comma-separated values?

Data frame A, which has 2 variables (mouse gene name and Ensemble ID):

head(dfA)
# gene EnsembleID
# 1 Hoxa9 ENSMUSG00000038227
# 2 Zfp128 ENSMUSG00000060397
# 3 Zfp853 ENSMUSG00000093910
# 4 Nr1h2 ENSMUSG00000060601
# 5 Nr1h3 ENSMUSG00000002108
# 6 Nr1h4 ENSMUSG00000047638

Data frame B has 3 variables (orthogroup, mouse gene Ensemble ID, and zebrafish gene Ensemble ID). Noted that in the picture below the column names are Mus_musculus.GRCm39.pep.all and Danio_rerio.GRCz11.pep.all, respectively.

![enter image description here][1]

I want to add a new column for "zebrafish Ensemble ID" into A as referenced information in B. So, the final data I need looks like this:
![enter image description here][2]

The point is, in data frame B, both Ensemble IDs' columns contain comma-separated values. This file is originally from the outputs of [orthofinder][3].

If there is only a single value in each row, I understand that we can easily combine those data frames using `merge` function in R. The length of data frame A is 1700 (= 1700 genes), so I would like to use some codes instead of search©&paste manually. Any suggestions or ideas would be really appreciated!

[1]: /media/images/bce57ec0-9892-475d-b704-e24fba59
[2]: /media/images/8c69cefd-2ae1-4c84-9ae7-a0104d18
[3]: https://github.com/davidemms/OrthoFinder

Answer: For dataframe manipulation, in general, you should look into the `dplyr` and `tidyr` packages, they offer endless possibilities if you learn to manipulate them (lots of practice will help). A good and concise cheatsheet is available [here][1].

Regarding this problem in particular, something like this should work:

library(dplyr)
library(tidyr)

dfA <- data.frame(gene = c("a", "b", "c", "d", "e"),
mus_id = c("mus001", "mus002", "mus003", "mus004", "mus005"))

dfB <- data.frame(mus_id = c("mus001", "mus002", "mus003,mus002", "mus004"),
zeb_id = c("dar001", "dar002,dar004", "dar003,dar002", "dar004"))

dfC <- dfB %>%
tidyr::separate_rows(zeb_id, sep = ",") %>%
tidyr::separate_rows(mus_id, sep = ",") %>%
dplyr::group_by(mus_id) %>%
dplyr::summarise(zeb_id = paste0(sort(unique(na.omit(zeb_id))), collapse = ','))

dfD <- dfA %>%
dplyr::left_join(dfC)

[1]: https://www.rstudio.com/wp-content/uploads/2015/02/data-wrangling-cheatsheet.pdf
Question: Hi,
This is a statistics question! I hope you may help me !
On a DNA sequence region, there are several mutation hits. I want to know if the hits distribution follow or not a random distribution.
In other words, imagine a target dart. How can get a p-value, to know if darts has been throw randomly or not .

![Exemple][1]

[1]: https://framapic.org/2fsQ36tnHGKL/LU4FryMZ2IDU.png

Answer: Both distributions are random, with the number of variants in a region as a discrete random variable. The opposite of random is deterministic, which would mean that the number of events is fixed, which is obviously not what you mean. Your second example shows a **uniform** distribution of events, while the first case should look like a deviation from uniform distribution, such that you would look for a test for deviation from uniformity. See here: https://stats.stackexchange.com/questions/25827/how-does-one-measure-the-non-uniformity-of-a-distribution
The simplest test proposed here is a Chi-square test.
Question: Hello,

After running [cellranger multi](https://support.10xgenomics.com/single-cell-gene-expression/software/pipelines/latest/using/multi) with single-cell VDJ-T sequencing data, I get the `cellranger` output for each sample. Is there a package that can do downstream analysis of single-cell V(D)J sequencing, like `scanpy` or `seurat` for single-cell RNA-seq?

Thanks a lot.

Answer: One option is Immunarch
There is also a dedicated immune repertoire sequencing community forum at that may be a good resource for you as well
Question: Hi everybody,

I am working on array CGH data from a set of tumor samples with potentially high normal contamination. I have performed the necessary pre-processing (background subtraction, within array normalization) on the tumor data, and run the circular binary segmentation on the resulting log ratio values.

I have seen that people often threshold the mean segment values following this step to call gains and losses (a typical value I have seen in literature is a minimum absolute value of 0.3 for log2ratio). However, this process assumes a minimum value of tumor purity, which may not be true for a subset of my samples.

I was wondering if anyone can suggest a more systematic way of approaching this problem in presence of high level (50-80%) of normal tissue contamination in the tumor sample.

Thank you!

Answer: In absence of purity estimates, thresholding is pretty much your only option. If your data is noisy then you might opt to focus on high-level amplifications and homozygous deletions in your report. Given another source of data like SNV calls, you can also use the allele frequencies detect loss of heterozygosity, which can support copy number calls for hemizygous losses at least.

Given purity estimates, you can rescale the log2 values with math, and then use thresholds more confidently.

For example, if your segmented aCGH calls are in SEG format (the output of DNAcopy), then you can use CNVkit to adjust the log2 values:

1. Load a sample into CNVkit's own format:

```
cnvkit.py import-seg Sample.seg -o Sample.cns
```
2. Rescale the log2 values using a purity estimate and optionally the ploidy:

```
cnvkit.py rescale Sample.cns --purity 0.45 --ploidy 2 -o Sample-rescaled.cns
```
3. Optionally, perform the thresholding by assuming 100% purity (now valid) and rounding log2 values to the nearest absolute integer copy number:

```
cnvkit.py call -m clonal Sample-rescaled.cns -o Sample-called.cns
```
Or by using hard cutoffs (in log2 scale):

```
cnvkit.py call -m threshold -t=-1.1,-0.4,0.3,0.7 Sample-rescaled.cns -o Sample-called.cns
```
4. Export the adjusted segments back to SEG:

```
cnvkit.py export seg Sample-called.cns -o Sample-called.seg
```
Or to another format:

```
cnvkit.py export bed Sample-called.cns -o Sample-called.bed
cnvkit.py export vcf Sample-called.cns -o Sample-called.vcf
```
Question: I have a vcf file with SNPs called using the Torrent Variant Caller and I want to calculate their Variant Allele Frequencies.

The issue is that there is both a DP (depth) and a FPD (flow depth) fields as well as an AO (Alternate Allele Observations) and a FAO (Flow Alternate Allele Observations) fields and the two often disagree by a lot.

Which pair should I trust? DP & AO or FDP & FAO?

Or should I just ignore the cases where they disagree?

Thanks

Answer: The fields FPD and FAO are after correcting for flow order of nucleotides, and so should be more accurate.

From https://ionreporter.thermofisher.com/ionreporter/help/GUID-A6D3AE1A-84C6-4B66-AA0D-69244576D2FE.html it looks like they don't use FPD in their allele freq calculations:

"AF=FAO/(FAO+FRO) and not FAO/FDP. This is because FDP may include reads that don't fit the flow space profile of any hypothesis; in such cases, FDP>=FAO+FRO and this is not used in allele frequency calculation"
Question: hi everyone,

I have just started learning genomics as a part of my bioinformatics degree and I've been introduced to linux for handling fastq files and using fast qc. can anybody suggest some good learning resources which is more inclined towards linux for genomics.

Answer: As you were told already, it may be a good idea to become Linux-proficient first. Many summaries of main Linux commands are available in 1-4 page formats if you Google `linux cheat sheet` or `linux command reference`.

Online courses, in lecture format:

- http://www.doc.ic.ac.uk/~wjk/UnixIntro/
- http://tldp.org/LDP/intro-linux/html/
- http://www.ee.surrey.ac.uk/Teaching/Unix/
Question: Hi there,

Fairly new to this area so will try to explain as clearly as possible. I have two RNA-Seq data sets - one corresponds to a series of cancer cell lines, the other to cell lines we are using as a model of 'normal' epithelium. The expression units of one data set are in FPKM the other RPKM. I think I superficially understand the difference between these units, in that one is used in mapping transcripts in single end sequencing the other paired end sequencing.

My question is are these units directly comparable? The analysis I wish to carry out is fairly straight forward - I have a predefined list of around 500 genes, and simply want to compare differences in expression between the non-cancer/cancer background, in terms of which of these genes they are expressing at all as well as the relative expression levels. For which transcripts are expressed I had intended to use any value over 0 (FPKM or RPKM) as denoting expression of a transcript, but am unsure if I can compare relative abundance.

I should add that I only have access to the raw data of one dataset, the other is as a results table sent by collaborators.

Answer: The difference between RPKM and FPKM is generally going to be much less than the difference between RPKM and RPKM for two different experiments with different methodologies. That is to say, since one is in RPKM and the other FPKM, the processing was probably different and thus they are not safe to compare regardless of the units. I recommend you try to get a hold of the original data and do the processing identically if you wish to compare results. Even then there may be other factors out of your control (wet-lab protocols) that make them incomparable.

You cannot convert between RPKM and FPKM, but in most cases they are extremely similar so if that were the only difference I doubt it would be significant.
Question: I am new to HISAT2 and trying to play around with its options. I found the following issues and would appreciate your advice:

1. I downloaded the genome_tran (for grch38) folder, and wanted to check out the option: --tmo. According to the [manual][1], it should "Report only those alignments within known transcripts.". However, this option allows <10% of my reads to align, while not using this option more than 80% of the reads could get aligned.

So I went further to examine the reads that could not get aligned after using --tmo, as far as I can manually tell, many of them are from coding regions. What do I miss here - any suggestions?

2. I also applied --rna-strandness F and --rna-strandness R respectively (with --tmo option on). My reads should be F. However, F and R gave me exactly the same output, with only diffice at " XS:A:" - one is +, the other gives - for a same read.

Thank you very much

[1]: https://ccb.jhu.edu/software/hisat2/manual.shtml

Answer: About `--tmo`, the issue has been [asked][1] and [reported][2] before, but it seems unresolved. It seems it works better with HISAT2 2.0.4.

About the `--rna-strandness`, I think this is the described and intended output:

> With this option being used, every read alignment will have an XS
> attribute tag: '+' means a read belongs to a transcript on '+' strand
> of genome. '-' means a read belongs to a transcript on '-' strand of
> genome.

Probably this information is used for downstream analyses, such as Stringtie and / or Ballgown.

[1]: https://www.biostars.org/p/305072/
[2]: https://github.com/infphilo/hisat2/issues/159
Question: Hi guys! I need help in setting a `for` loop in R (I'm quite new programming in R). I would like to add the same header to all the files that match a concrete pattern inside a folder.

To get the list of files I'm using the following code:

filelist <- list.files(pattern = "DESeq2_result*")

And this is the `for` loop I am trying to implement:

for (i in seq_along(filelist)) {
names[[i]] <- a
out [i]
}

where `a` is a vector that I defined with the names of the different columns:

a <- c("gene_id", "baseMean", "log2FC",
"SD", "WaldStatistic", "pval", "padj")

If you have any tutorial/page to help me to learn and practise my ability to code functions and loops in R would be so much appreciated.

Thank you so much in advance!

Jordi

Answer: As I understood, you want to open some files containing the pattern "DESeq2_result*" as data.frames in a list, assigning new names to the columns of these dataframes based on the `a` vector, and maybe name the elements of the list based on the files name? You can do it that way:

filelist <- list.files(pattern = "DESeq2_result")

a <- c("gene_id", "baseMean", "log2FC",
"SD", "WaldStatistic", "pval", "padj")

data<-list()

for (i in 1:length(filelist)) {
data[[i]]<-read.table(filelist[i])
colnames(data[[i]])<-a
}

names(data)<-filelist

for (i in 1:length(data)) {
output_name <- paste0("colnames_",names(data[i]))
write.table(data[[i]], output_name, quote = F, row.names = F, col.names = T, sep="\t")
}

EDIT: added the write.table part to take your comment into account in order to export the imported files with the column names.
Question: Hello,

Can anyone suggest a simple R package to perform GO term analysis. Specially I am looking for a package that can produces publication quality figures and heatmaps regarding GO term results.

Thanks in advance

Answer: Check bioc package "clusterprofiler".
https://guangchuangyu.github.io/2016/01/go-analysis-using-clusterprofiler/

For heatmaps of GO enrichment, check "revigo". They provide the output as R script which you can load in R and manipulate colors etc. to suit your need.
Question: Hello, I have a question for you guys. The perl script rfam_scan.pl uses the program cmsearch to look for a CM model. Wouldn't be better to use cmscan instead of cmsearch? According to the Infernal manual, we should use cmscan if we have a set of sequences and we wanna know the corresponding or closest CM. Don't you think so?

Answer: As far as I recall, cmscan (like hmmscan) and cmsearch (like hmmsearch) return the exact same results. The big difference is that the former can't be parallelized without special effort.

Edit. http://selab.janelia.org/people/eddys/blog/?p=424
Question: I have a list of SNPs and their exact locations :

SNP_nr Chromosome Location
1 1 11088169

How can I get their associated *rsid* programmatically or otherwise?

Answer: many answers on biostars , search for UCSC/mysql/api.

eg:

Is there such a thing as a UCSC API?

direct acces to mysql: (slow, because it doesn't use the 'bin' index)

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -Dhg19 -e "select name from snp137 where chrom='chr1' and chromStart=10144"
+-------------+
| name |
+-------------+
| rs144773400 |
+-------------+
Question: Hi all,

I have a scRNA-seq dataset, which has 6 patients and each patient has 2 sample types (normal, tumor). So I have **12 folders** and each folder contain its scRNA data: `barcodes.tsv.zip`, `features.tsv.zip`, `matirx.mtx.tsv.zip`. I would like to use this data to practice and learn the `Seurat (version 4.0)` r package workflow.

I would like to process all the 12 samples into cluster and downstream analysis. For the preprocessing, I now have no idea which point I should merge all the datasets.

Should I **read and QC each folder** (e.g patient1-normal -> patient1-tumor -> patient2-normal .... so on), and then merge all data? I am still learning the scRNA seq, I believe each data will have different number of genes and cells, I need to normalize it after merging all the data.

Thank you for your help and advice.

Answer: Yes, you should QC each sample individually (you can automate this of course but be sure to manually inspect each QC plot). This is at least what I do. There are approaches like 3xMAF to automatically filter data, and I use them in support of what my eyes see in the QC plots, but I personally think that the human brain is an awesome diagnostic tool so just looking at the plots makes a lot of sense to see whether the automated cutoffs are actually meaningful. With skewed data they might not be. And then for clustering you should integrate your data. Please read the Seurat integration vignette. I suggest though that you download some of the Seurat example data to practice with because with 12 samples it will take both time and memory. Better to go easy with a smaller dataset or just two patients to start with. As an alternative to Seurat with lots of code to follow is https://osca.bioconductor.org/, I recommend it as documentation is (for me) much clearer. It has a section on multi-sample integration as well.
Question: Greetings,

I have genetic maps for all the chromosomes as follow:

Chromosome	Position(bp)	Rate(cM/Mb)	Map(cM)

chr1	55550	2.981822	0.000000

chr1	82571	2.082414	0.080572

chr1	88169	2.081358	0.092229

...

And I want to convert Rate(cM/Mb) to a probability value (the recombination probability), does anyone have an idea about how to do that?

Thanks,

Answer: Here is the answer to my question :)

I need to take last column "Map(cM)" and by using the last formula in https://en.wikipedia.org/wiki/Centimorgan I can convert it to recombination probability.
Question: Hi,

I've got an error using Picard MarkDuplicates function. It worked fine on my data before I add an other step in my pipeline, where I use samtools view -F4 -h to remove unmapped reads in the alignement. If I do this, an error occured during Picard treatment :

> Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 6, Read name FCC63KDACXX:3:1115:14357:43781#GCGGAACT_TCTTTCCC, Mate unmapped flag should not be set for unpaired read.

Do you have any ideas of how to fix that ? I understand the error but I don't know how to make it works. Previously I 'manually' removed unmapped reads (reads for which an * appears in the 3rd column in the sam files) but I read that's not the good way to do that and that the best way is to use samtools with the command line above.

Thanks.

Answer: Doing -F 4 will result in reads that were not mapped being removed from the file - but if their mate was mapped, the mate will be left in there as a *singleton*. Picard expects these singletons to look like singletons - not to look like half a pair, with the other half just mysteriously gone.

I imagine that running [FixMateInformation][1] will sort this out, otherwise you can use -F 12 instead of -F 4 when filtering the BAM. I can't remember if samtools does AND or OR when using bitflags, so that will either remove both reads in a pair if the read AND the mate is unmapped, or both reads in a pair if the read OR the mate is unmapped (what you probably want).

[1]: http://broadinstitute.github.io/picard/command-line-overview.html#FixMateInformation