cannin/biostars_qa huggingface dataset:
Question: I have a set of RNAseq data from 4 wt, 10 hets and 2 ko. I am trying to see if more replicates in hets leads to more degs in het vs ko than wt vs ko. Hence, I am here trying to take only 4 random replicates of hets from the DeseqDataSet and run analysis multiple times.

In each run, only difference is which hets are in the DESeqDataSet. All runs have have the same wt and het. Even then, my DEGs for wt vs ko are different, in each run. How is that?

random_het <- function(dds_lcr){
minusHet <- sampleTable %>% filter(genotype %in% "KO.WT") %>% select(embryo) %>% distinct() %>% slice_sample(n =6) %>% as.list %>% unlist %>% paste0(., collapse = "|")
dds_4het <- dds_lcr[, -grep(minusHet, colnames(dds_lcr))]
dde_4het <- DESeq(dds_4het)
WtVsKo4het <- results(dde_4het, name = "genotype_WT_vs_KO", lfcThreshold = 0, alpha = 0.05)
HetVsKo4het <- results(dde_4het, name = "genotype_KO.WT_vs_KO", lfcThreshold = 0, alpha = 0.05)
sum.data <- map(list(WtVsKo4het, HetVsKo4het), getSig, padj = 0.05) %>% map(., nrow) %>% as.data.frame %>% set_names(c("WtVsKo4het", "HetVsKo4het"))
names(sum.data) <- c("WtVsKo4het", "HetVsKo4het")
data <- list(sum.data, embyos = colData(dds_4het))
return(data)

}

randomise_plot <- function(n){
randomiseHet <- purrr::rerun(n, random_het(dds_lcr))
df_dge <- sapply(randomiseHet, "[[", 1) %>% t %>% as.data.frame
plot(df_dge)
p1 <- recordPlot()
dds_names <- sapply(randomiseHet, "[[", 2)
return(list(p1, dds_names))
}

randomise_plot(n=10)

The plot shows in two runs there are different number of DEGs for both comparison. [plot][1]

[1]: https://drive.google.com/file/d/1yVWnPFNtWsb3r_hJXqmLqGCfdst-Pifb/view?usp=sharing

Answer: Even if you compare only WT and KO for differential expression, DESeq2 internally integrates all samples (with the included hets in your case) to estimate dispersion. As dispersion is a measure of within group variability, changing what sample are included in a group (hets) will change the dispersion estimate, which will change the confidence of calling a gene differentially expressed.

To give you an example, consider the two following situations:

# expression for gene x (CASE A)
wt_1 100
wt_2 110
KO_1 200
KO_2 205
het_1 100
het_2 1000

# expression for gene x (CASE B)
wt_1 100
wt_2 110
KO_1 200
KO_2 205
het_3 200
het_4 200

In the first case, the dispersion estimate will be very high because the within group variability is huge in the het group. Therefore, the gene will not be significantly differentially expressed, even between condition wt and ko (within groups variability >>> between groups variability). In the second case, it most likely will be differentially expressed, because the within group variability and the dispersion estimate are much smaller.

For those reasons, it is best to include all samples in the differential expression analysis (unless there is a clear outlier/artifact). This usually provides the best dispersion estimate.
Question: **Question**:

What is the simplest/fastest/safest way to rename a few INFO tags in a VCF file (especially if the file is already bgzipped)?

**Background**:

I'm trying to merge multiple VCF files, but some of them have conflicting INFO tag names (e.g. AC, AN, AF, etc.). I want to rename some of the INFO tags in the original (unmerged) files so they don't conflict with each other once they're merged.

**What I've tried**:

BCFtools is nice in that it's safe, flexible and works very quickly with bgzipped VCF files. I figured this simple task would be builtin, but I couldn't find it anywhere. I tried using `bcftools reheader` to rename some tags, but this only modifies the VCF header. The tags in the data remain the same. I'm open to using other tools (VCFtools, awk, sed, etc.), but I would prefer a method that is compatible with bgzipped files (if possible).

Answer: Here's the best I could come up with for now:

```
bcftools view -O v variants.vcf.gz \
| sed -e 's/\([;=[:space:]]\)EA_AC\([,;=[:space:]]\)/\1EVS_EA_AC\2/' \
-e 's/\([;=[:space:]]\)AA_AC\([,;=[:space:]]\)/\1EVS_AA_AC\2/' \
| bcftools view -O z -o variants.re-tagged.vcf.gz
```

**Explanation**:

1. I'm using `bcftools view -O v` to open my compressed VCF file (in this case `variants.vcf.gz`). This command is nice because it also supports opening uncompressed VCF files (e.g. `variants.vcf`) as well as compressed/uncompressed BCF files (e.g. `variants.bcf.gz` or `variants.bcf`, respectively) without having to specify any extra flags.
2. Then I'm using `sed` to rename my INFO tags (both in the header and the data). In this case I'm renaming `EA_AC` and `AA_AC` to `EVS_EA_AC` and `EVS_AA_AC`, respectively.
3. Finally I'm using `bcftools view -O z` to write my output to a compressed VCF file (i.e. `variants.re-tagged.vcf.gz`). This command also allows you to write to an uncompressed VCF or compressed/uncompressed BCF. If you leave this step off then the uncompressed output will be written to standard out.
Question: Dear all,

although this question has been asked before (~ 2 years ago), if I may ask it again, as I am looking for some updated workflows, strategies, ideas :

"what would be the acceptable pipeline for processing WES or WGS data from primary tumors that do not have matched germlines sequences ?"

many thanks,

-- bogdan

Answer: MuTect and FreeBayes for SNVs are probably the most used tools for tumor-only analyses. Copy number, if you have some normals, you can try CNVkit, cn.mops (I think) and GATK4beta. If you need purity/ploidy and want to classify variants as germline vs somatic, there isn't much out yet, so we wrote our own tool PureCN (Bioconductor, try the much improved developer version if you like).
Question: Good evening,

I followed the tuxedo protocol as described [here][1]. With the "normal" outputs that cufflinks provide (I'm searching on cuffdiff output files, but tell me if I'm wrong), how could I determine which of isoforms detected are actually novel ones?

I'm pretty new on this field, and I've found in many posts that using Cufflinks is the answer to find novel isoforms, or differential splicing too. Now that I followed up the "instructions" I'm stuck on trying to interpret the results, after searching all over I could. Searching on cufflinks manual didn't help me either. :/

Can someone explain me or provide me some tutorials/papers/etc, that might help me on this topic please?

Thanks in advance.

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC3334321/pdf/nihms-366741.pdf

Answer: Here are two tutorials on RNA-seq, both that involve discussion of Cufflinks that you might find generally useful:

- http://www.rnaseq.wiki
- http://rnaseq.uoregon.edu/

With regard to your specific question. You might find the Cuffcompare and Cuffmerge commands and their documentation useful.

- http://cole-trapnell-lab.github.io/cufflinks/cuffmerge/index.html
- http://cole-trapnell-lab.github.io/cufflinks/cuffcompare/index.html

Cuffcompare has this description:

> After assembling a transcriptome from one or more samples, you'll probably want to compare your assembly to known transcripts. Even if there is no "reference" transcriptome for the organism you're studying, you may want to compare the transcriptomes assembled from different RNA-Seq libraries. Cuffcompare helps you perform these comparisons and assess the quality of your assembly.

You can run cuffcompare directly but in recent versions of cufflinks, the cuffmerge command will do this for you. You should have '.tracking' output files. Refer to the 'class code' column of this file. The definition of "transfrag class codes" is provided at the bottom of the Cuffcompare manual page. For example, class "j" is described as "Potentially novel isoform (fragment): at least one splice junction is shared with a reference transcript"
Question: Hi, I've noticed that the Wilcoxon rank-sum test is used by many R packages to calculate differential expression. Is this a good test to estimate that when you compare it to other packages like DESeq2 and edgeR which calculate differential expression with many more factors involved and would that not be better than just a Wilcoxon test?

Answer: I agree with Devon here except I'd put the number even lower than 100 for trusting Wilcoxon. Rank based tests will find differences with sufficient power once you have, say dozens of cells with non-zero counts in each group that you want to compare.

We used Wilcoxon tests with Gibbs samples to also account for uncertainty from multimapping reads in Swish:

https://academic.oup.com/nar/article/47/18/e105/5542870

Jump to "scRNA-seq simulation" and "scRNA-seq of mouse neurons" for an examination of how Wilcoxon can be used effectively for comparing moderate-sized groups of cells (e.g. ~20 or higher) even in the presence of uncertainty from multi-mapping reads.
Question: Hi.

I generated a heatmap with several clusters using "deepTools" something like below, and I want to know the names of elements in each cluster. But, it seemed no tool existed among the "deepTools suit".

Can anyone help me?

Thank you!

![enter image description here][1]

[1]: http://deeptools.readthedocs.io/en/latest/_images/ExampleHeatmap3.png

Answer: If you use the `--outFileSortedRegions something.bed` option then a BED file will be produced with a column noting the cluster that each gene/region/feature is in.
Question: Hi all,

I want to normalize the Illumina Infinium methylation microarray data by lima, or by minfi or other programs you would recommend.

The online available GEO... has: *family.soft data, matrix file, and RAW data like: v1.1.bpm, v1.2.bpm, *.csv, manifest/header_descripton.xls. and RAW data with:
TargetID	Proband393_mother.Signal_A	Proband393_mother.Signal_B	Proband393_mother.Detection Pval	Proband393_father.Signal_A	Proband393_father.Signal_B	Proband393_father.Detection Pval	and so on....

I checked the lumi tutorial and the package required IDAT data.

Could anyone tell me what data should I use to the normalization? I want to get the Methylation fractions (of Beta value) and the genomic locations of them.

I would be very grateful for help.
I am new in epigenetics and still learning.

Dorota

Answer: Dear Dorota,

The normalised values for your project of interest, GSE89353, can be obtained by running this code:

library(Biobase)
library(GEOquery)

gset <- getGEO("GSE89353", GSEMatrix =TRUE, getGPL=FALSE)
if (length(gset) > 1) idx <- grep("GPL13534", attr(gset, "names")) else idx <- 1
gset <- gset[[idx]]

boxplot(exprs(gset)[,1:50])

--------------------------------------
-------------------------

You can then do further tests with this data using standard R functions, e.g., difference in mean between different groups, or use other packages like:

- https://bioconductor.org/packages/release/bioc/vignettes/minfi/inst/doc/minfi.html
- missMethyl

Kevin
Question: I often see the word depth in the manuals of the tools for NGS, what is its meaning ?

thanks.

Answer: Eric gives the correct answer for depth (of coverage). I think confusion in this area stems not from the term "depth" but from the term "coverage". Coverage now appears to have 3 meanings:

the theoretical "fold-coverage"
of a shotgun sequencing experiment:
number of reads * read length /
target size

the theoretical or
empirical "breadth-of-coverage" of
an assembly: assembly size / target
size

the empirical average
"depth-of-coverage" of an assembly:
number of reads * read length /
assembly size

(1) and (3) are not the same because of sequencing error & unclonable/unmappable regions of the genome. Lander-Waterman theory deals with the relationship between (1) and (2).
Question: Hi,

I have a data set which includes A, T,C,G.
When I use read.table command for this data set, I see TRUE's in terms of T alleles. It happens if the entire column just includes T. How can I fix that problem?

Thanks.

Answer: `read.table()` tries to guess the class of the input data and will sometimes be mislead.

> read.table(text=("A C G T"))
V1 V2 V3 V4
1 A C G TRUE

> summary(read.table(text=("A C G T")))
V1 V2 V3 V4
A:1 C:1 G:1 Mode:logical
TRUE:1
NA's:0

It is possible to specify in advance the class of the columns; see `?read.table` for details.

> read.table(text=("A C G T"), colClasses = "character")
V1 V2 V3 V4
1 A C G T

Related to Chris's answer in this particular case, `stringAsFactors` will not solve the problem.

> read.table(text=("A C G T"), stringsAsFactors = FALSE)
V1 V2 V3 V4
1 A C G TRUE

Note that there are other cases where `T` may be coerced to `TRUE` instead of `"T"`. In particular, pay attention that [there is one gene whose symbol is _T_](https://en.wikipedia.org/wiki/Brachyury) !
Question: Hi all,

Let's imagine I'm developing a Python tool that can predict clinical phenotype for the gene YFG (Your Favorite Gene) which is located in chr1:1000-2000. It does so by first creating a VCF file from an aligned BAM file using `bcftools` (more specifically via the Python package `pysam` -- i.e. `pysam.bcftools`) and then somehow mapping variant records to a known phenotype.

For variant calling the tool currently requires two input files from a user: a BAM file and a reference FASTA file. However, I want to update the tool so that it already includes a "mini" FASTA file that only contains the sequence of YFG which is 1 kb long. This means users are now required to provide their BAM only. I can't just put an entire reference FASTA file to the tool because it is usually huge (~3 GB) so the tool won't be portable anymore (portability matters to me a lot). Including mini FASTA will also improve the reproducibility of said tool because currently different users have slightly different FASTA files (e.g. different contig names such as `chr1` vs. `1`).

That being said, my question is: Is this feasible? I keep thinking there got to be a simple way to achieve my goal, but I can't think of a good way to actually pull this off.

For instance, if I were to simply use mini FASTA as described above, `bcftools` will return most positions in chr1:1000-2000 as variants because every position has been left shifted by 1000. I know this because I have already tried it :(

I still think above approach could work if, for example, there is a way to tell `bcftools` that the sequence in mini FASTA starts at chr1:1000. But no luck thus far. Any advice and recommendation would be highly appreciated.

Another possible workaround might be to re-align sequence reads from BAM to mini FASTA and then perform variant calling with mini FASTA. I actually think this will work quite well, but the major problem is that there is no Python package that can perform efficient and accurate read alignment (e.g. `bwa-mem`). Of course, I can have users install a read aligner separately, but if I can I really want to package everything within the tool.

Thanks for reading this lengthy post!

Answer: In normal operation, the coordinates of the BAM and FASTA file need to match.

One workaround could be to extract a subset of the BAM file, then post-process that BAM file to shift the POS column to match the mini-FASTA (subtract the offset from each POS field - other corrections may also be needed) then use the mini BAM file with the mini FASTA
Question: Hello All,

I'm new to the forum and I'm looking for a bit of help. I'd like to pull down the Refseq fasta protein sequences of a gene (all isoforms), however I'd like to do that with several species. For example, I'd like to get the gene ACTB for human (which has only 1 isoform NP_001092.1), but also get sequences for mouse (NP_031419.1), and macaca fascicularis ( NP_001271954.1). This is easy for an example since each only has one isoform, and the monkey assembly is well annotated for this gene so it doesn't have the XP_ type sequence it usually does.

This procedure is a bit tricky because sometimes the gene I'm pulling may not exist in one of the other species (i.e. not present in mouse), and importantly I need them all to be from the same reference source build (i.e. should all be RefSeq or all Ensembl). I have a very strong preference for RefSeq though.

Here's what I've tried using other methods in R, but no success:

library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(rentrez)
symbols <- c("^ACTB$")
grep(symbols, keys(x = org.Hs.eg.db, keytype = "SYMBOL"), ignore.case = T, perl = T) #test if gene found in human
grep(symbols, keys(x = org.Mm.eg.db, keytype = "SYMBOL"), ignore.case = T, perl = T) #test if gene found in mouse
keys.h <- keys(x = org.Hs.eg.db, keytype = "SYMBOL")[grep(symbols, keys(x = org.Hs.eg.db,
keytype = "SYMBOL"), ignore.case = T, perl = T)]
keys.m <- keys(x = org.Mm.eg.db, keytype = "SYMBOL")[grep(symbols, keys(x = org.Mm.eg.db,
keytype = "SYMBOL"), ignore.case = T, perl = T)]
#map HGNC to ID#
id.h <- mapIds(x = org.Hs.eg.db, keys = keys.h, keytype = 'SYMBOL', column = 'ENTREZID')
id.m <- mapIds(x = org.Mm.eg.db, keys = keys.m, keytype = 'SYMBOL', column = 'ENTREZID')
select(x = org.Hs.eg.db, keys = id.h, columns = c("ENTREZID", "REFSEQ", "SYMBOL"))
select(x = org.Mm.eg.db, keys = id.m, columns = c("ENTREZID", "REFSEQ", "SYMBOL"))
#pull Human protein data
id.link <- entrez_link(dbfrom = "gene", id = id.h, db = "protein")
id.prot <- id.link$links$gene_protein_refseq #this is where no protein is found
id.seqs <- entrez_fetch(db = "protein", id = id.prot, rettype = "fasta")

As you can see I'm able to retrieve human and mouse decently, but unfortunately I can't figure out a way to pull the monkey protein sequences. I can't find an equivalent "org.Mf.eg.db"database for the monkey, that would pull the NP_001271954.1 sequence. Can anyone help me figure this one out?

Thank you in advance!

Answer: well for any poor sod that comes across this post and wanted to know, I solved my own problem by utilizing the BSGenomes package and GenomicFeatures package:

> library(GenomicFeatures) cyno.5.0 <- makeTxDbFromUCSC(genome =
> "macFas5", tablename = "refGene") saveDb(x = cyno.5.0, file =
> "TxDb_cyno_UCSC_5_refGene")
This will build a TxDb object that can be used to get the sequences for translation. A bit convoluted but works for the less used organisms. This could be amended for another species by simply scanning other available builds (uses rtracklayer package loaded as dependency):

ucscGenomes() #scroll through to find relevant organisms
Question: In [this paper][1] , section "Affymetrix platform", the below paragraph has been written:

*Affymetrix was the first to commercially produce SNP arrays, nearly a decade ago. The HuSNP assay, initially prototyped in Wang et al. (15), was designed to genotype 1494 SNPs on one chip. Subsequent versions increased stepwise from 10 000 to 100 000 to 500 000, and finally to nearly one million SNPs in the current release. Every SNP site is interrogated by a set of probes that are each 25-nt long. A probe is designed to be complementary, or very nearly complementary, to a portion of the DNA sequence harboring the SNP site (Figure 2a). In the first few versions of the array, each SNP was interrogated by between 24 and 40 distinct probe sequences, forming a probe set. Within a set, each probe is associated with either allele A or allele B. Additionally, each probe is either a perfect match (PM; perfectly complementary to one of the target alleles), or a mismatch (MM; identical to a perfect match probe except that the center base is altered so as to be perfectly complementary to neither allele). The idea of the mismatch probe comes from mRNA expression arrays (16), and their purpose is to measure background noise. The scheme yields quartets comprised of four types of probes: PMA, MMA, PMB and MMB. The computational goal is to convert these 8–10 probe quartet intensity measures from raw array data into a genotype inference—AA, AB or BB.*

In the last sentence, it has been written that "8–10 probe quartet". Why "8-10"?

[1]: http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2715261/

Answer: For each of the types (PMA, MMA, PMB, MMB) there are upto 10 probes. (hence the 24-40 distinct probes earlier in the paragraph)
The earlier SNP arrays also had shifting positions around the SNP
if p is the basepair position of the SNP, then there were probes between p-12 and p+12, but also p-11 to p+13 and p-13 to p+11 (all 25 nucleotides long)

I don't have experience with the earliest type with 1500 SNPs, but IIRC the 10k SNP array had 5 probes for each type, each shifted by 1 position. ( p-14 to p+10, p-13 to p+11, p-12 to p+12, p-11 to p+13, p-10 to p+14)
Later Chiptypes all have fewer probes per SNP, and at a certain time the MM probes were discarded alltogether
Question: Hi dudes,

I want to count the number of my reads that have been mapped to each gene, then I downloaded featureCounts from http://sourceforge.net/projects/subread/files/subread-1.4.6-p5/, untarred that then using **my accepted_hits.sam** and **genes.gtf** as input, I replaced the command from documentation like below:

**Summarize a single-end read dataset using 5 threads**:

featureCounts -T 5 -t exon -g gene_id -a annotation.gtf -o counts.txt mapping_results_SE.sam

```
[izadi@lbox161 Downloads]$ ls
subread-1.4.6-p5-source subread-1.4.6-p5-source.tar.gz
[izadi@lbox161 Downloads]$ cd subread-1.4.6-p5-source
[izadi@lbox161 subread-1.4.6-p5-source]$ ls
LICENSE README.txt annotation doc ouput src test
[izadi@lbox161 subread-1.4.6-p5-source]$ cd src
[izadi@lbox161 src]$ ln -s /usr/data/nfs6/izadi/angel/ouput/
[izadi@lbox161 src]$ featureCounts -T 5 -t exon -g gene_id -a genes.gtf -o counts.txt accepted_hits.sam
bash: featureCounts: command not found...
```

After untaring there no bin file there? do you know why?

[izadi@lbox161 izadi]$ tar -zxvf subread-1.4.6-p5-source.tar.gz

May you please read my command and detect my fault here?

Thank you

Answer: Enter the src subdirectory under the home directory of the package and then issue the following command to build it on a Linux/unix computer:

make -f Makefile.Linux

then the bin will be created and Enter that subdirectory under the home directory, now run your command

as simple as eating a piece of cake but not cheese one!
Question: Dear all,

lets assume we have reference genome and its mapped reads:

ref: ACGTCTTCGTCAATG
read1: ACGA
read2: ACGA
read3: ACGC
read4: ACGC

my variant caller should generate something like: **het A/C - 50%**. Dou you have any experiences with this calling? What caller can handle this problem? I found that some callers just skip this site.

Best,

Paul.

Answer: GATK's HaplotypeCaller would call this:

| CHROM | POS | REF | ALT | GT |
|-------|----------|-----|-----|------|
| chr | 1 | T | A,C | 1/2 |

So for GT the ref (T) is 0, the first alternate (A) is 1 and the second (C) is 2, so

TT = 0/0 = WT

AA = 1/1 = Homozygous variant first alternate

CC = 2/2 = Homozygous variant and the second alternate

TA = 0/1 = Heterozygote with REF and first alternate allele

TC = 0/2 = Heterozygote with REF and second alternate allele

and what you described above:

A/C = 1/2 = Heterozygote with both alternate alleles