cannin/biostars_qa huggingface dataset:
Question: Hello,

I'm trying to remove unwanted variables using RUVSeq and SVA but i'm meeting an issue with both libraries.
with RUVSeq I get this "list" error:

> for(k in 1:4) {
+ set_g <- RUVg(x = gg, cIdx = house_keeping_genes, k = k)
+ DESeq2::plotPCA(set_g, col=as.numeric(coldata$group))
+ }

```Error in Ycenter[, cIdx] : invalid subscript type 'list'```

while with SVA I get this:

> res <- svaseq(norm_by_count_per_million, mod, mod0)

```Number of significant surrogate variables is: 46```
```Iteration (out of 5 ):Error in density.default(x, adjust = adj) : 'x' contains missing values```

I've found a couple of references but haven't been able to fix these issues.
Full reproducible example bellow, any help would be much appreciated!

library(SummarizedExperiment)
library(RUVSeq)
library(sva)

# http://duffel.rail.bio/recount/v2/TCGA/rse_gene_bladder.Rdata
load('data/rse_gene_bladder.Rdata')

unique(rse_gene$gdc_cases.diagnoses.tumor_stage)

colData(rse_gene)$group <- NA
rse_gene$gdc_cases.diagnoses.tumor_stage == "stage iv"
colData(rse_gene)$group[rse_gene$gdc_cases.diagnoses.tumor_stage == "stage i"] <- "early"
colData(rse_gene)$group[rse_gene$gdc_cases.diagnoses.tumor_stage == "stage ii"] <- "early"
colData(rse_gene)$group[rse_gene$gdc_cases.diagnoses.tumor_stage == "stage iii"] <- "late"
colData(rse_gene)$group[rse_gene$gdc_cases.diagnoses.tumor_stage == "stage iii"] <- "late"

keep <- !is.na(rse_gene$group)
rse_gene <- rse_gene[, keep]
table(rse_gene$group)

counts <- assay(rse_gene, "counts")
norm_by_count_per_million <- sweep(counts, 2, FUN="/", colSums(counts)) * 10^6

coldata <- colData(rse_gene)

HK_genes <- read.table(url("https://m.tau.ac.il/~elieis/HKG/HK_genes.txt"))
house_keeping_genes <- intersect(rowRanges(rse_gene)$symbol, HK_genes$V1)
gg <- ceiling(norm_by_count_per_million)
rownames(gg) <- as.vector(rowRanges(rse_gene)$symbol)

par(mfrow = c(2, 2))
for(k in 1:4) {
set_g <- RUVg(x = gg, cIdx = house_keeping_genes, k = k)
DESeq2::plotPCA(set_g, col=as.numeric(coldata$group) )
}

mod <- model.matrix(~group, data=coldata)
mod0 <- model.matrix(~1, data=coldata)

res <- svaseq(norm_by_count_per_million, mod, mod0)

Answer: The `?RUVg` help:
> cIdx
A character, logical, or numeric vector indicating the subset of genes to be used as negative controls in the estimation of the factors of unwanted variation.

But `class(house_keeping_genes)` returns `list`, hence `as.character(house_keeping_genes)` should work.

I would recommend normalizing data with DESeq2 or edgeR, not with this plain per-million scaling. Here is a full DESeq2+RUVseq analysis from the DESeq2 author: https://github.com/mikelove/preNivolumabOnNivolumab/blob/main/preNivolumabOnNivolumab.knit.md
Question: I have many VCFs from different samples. If I merge these into a single VCF using `vcftools` (`vcf-merge`), the samples where a variant wasn't called are labeled as missing that variant. Instead, I want the VCF to show that the sample has the reference allele (safe to assume in my application).

Is there a way to call missing variants in a VCF as the reference allele? What tools can I use to do this?

**EDIT:**

The sequences were originally variant called using FreeBayes (through the LongRanger pipeline).

**RE-EDIT:**

Turns out I can simply use the `--ref-for-missing` flag in `vcf-merge` to achieve this. Problem solved.

**RE-RE-EDIT:**

Using `--ref-for-missing` flag in `vcf-merge` does of course not give the variants any annotation, like depth and genotype quality.

Answer: ( The best way is to call all the BAMs in the same command, to get a multi-sample VCF)

I've written two tools related to your question:

* **VcfNoCallToHomRef** http://lindenb.github.io/jvarkit/VcfNoCallToHomRef.html will arbitrarily convert the `./.` to `0/0` .

* **FixVcfMissingGenotypes** http://lindenb.github.io/jvarkit/FixVcfMissingGenotypes.html : will re-use the BAM to test the depth at each NO-CALL site. It's slow;
Question: Hi,

I heard about a tool which lets you search for papers based on a set of genes that you input and it returns papers that have mentioned a significant set of those genes.

Does anyone know what tool this is?

Thanks

Answer: In R you can do this using Org.HS. Below is an example:

```
source("http://bioconductor.org/biocLite.R")
biocLite("org.Hs.eg.db")
library("org.Hs.eg.db")

# mapped_genes are all the genes that org.HS.egPMID covers (the HS refers to homosapien). These genes are in Entrez format.

# entrez2Pmid is a list taking entrez gene ids to a vector of Pmids.

mapped_genes <- mappedkeys(org.Hs.egPMID)
entrez2Pmid <- as.list(org.Hs.egPMID[mapped_genes])

# now if you have an entrez gene id you can look up the relevant papers by doing

entrez2Pmid[myEntrezGene]
```

[Here][1] is a pastebin of the code since I stink at using the biostars code formatter

There is a tutorial [here][2].

Also you can just use http://idconverter.bioinfo.cnio.es/IDconverter.php. Although I've found that sometimes the results seem to be a subset of the results you get using the method outlined above.

[1]: http://pastebin.com/JCa6L7sB
[2]: http://www.bioconductor.org/packages/release/data/annotation/manuals/org.Hs.eg.db/man/org.Hs.eg.db.pdf
Question: Hi all,
I have a file like this

group1 group2 group3 group4
ax as we aw
as we rt ty
aw aq yu pl
aq qw oo se

I need to count unique (pairwise) one and save as a matrix form like this

group1 group2 group3 group4
group1 4 2 0 1
group2 2 4 1 0
group3 0 1 4 0
group4 1 0 0 4

I used R intersect function for this but I could not able to save the output as a matrix form.

Answer: How about double loop:

# example data
df1 <- read.table(text = "
group1 group2 group3 group4
ax as we aw
as we rt ty
aw aq yu pl
aq qw oo se
", header = TRUE, stringsAsFactors = FALSE)

# double loop: get length of intersect
sapply(df1, function(i)
sapply(df1, function(j) length(intersect(i, j))))

# group1 group2 group3 group4
# group1 4 2 0 1
# group2 2 4 1 0
# group3 0 1 4 0
# group4 1 0 0 4
Question: Hi,

I have a Python script which takes as input a list of SNPs, then calls tabix to retrieve .VCF file and does some stuff with it.
The problem is that i need to manually add (using UCSC genome browser web*) the coordinates of the SNP for tabix to work properly.

Does any one know any simple python API / function / module to easy add this information without having to perform the search manually?

Thanks!

Answer: You can do this with cruzdb

from cruzdb import Genome
import sys

fname = sys.argv[1]
hg19 = Genome('hg19')
snp135 = hg19.snp135

for rs in (l.split()[0].rstrip() for l in open(fname) if l.startswith("rs")):
print snp135.filter_by(name=rs).first()

which takes a file with the rs id's in the first columns and prints out a bed file with the chromosome locations from dbsnp 135.
Question: Hi, I'm very confused about the IGH gene. So on NCBI it looks to be a single IGH@ gene or locus. however

the region seems to cover many other IGH genes IGHD1, G1, exc...
this is how it looks like, [browser][1]

Why I'm interested is because of gene fusions where the 3' gene is written as IGH. However, this is ambiguous because this region seems to cover many variants. My questions are the following.
1. is IGH a single gene or locus containing many small variants. Or if this a gene that gets post prossessed post joining different segments?
1b. do these IGH "variant" contain exons/intron boundaries.
2. On the CCDS website I cannot seem to be be able to locate any IGH genes at all. What I want is to have the coordinates and sequences for each variant.
3. Is there a place I can get the the sequences for each segment and coordinate, best if I can find it on CCDS.

thanks in advance.

[1]: https://genome.ucsc.edu/cgi-bin/hgTracks?db=hg38&lastVirtModeType=default&lastVirtModeExtraState=&virtModeType=default&virtMode=0&nonVirtPosition=&position=chr14%3A105586437%2D106879844&hgsid=1475821023_Ps5LU7lpkag5rANMUtJpmWA2WtQW

Answer: 1.The IGH locus contains many V, D, and J gene segments, as well as constant region genes. The exact number of V, D, and J genes will depend on the species. For example, cows have 12 V, 16 D, and 4 J regions at the IGH locus. Canonical antibodies are composed of two heavy chains (HC) and two light chains (LC). B cells make antibodies. During B cell development, a V, a D, and a J gene combine to form the variable region of the HC. The variable region is joined to a constant region to make a full length HC. See the Wikipedia link below for a description of V(D)J recombination:

[V(D)J recombination][1]

1b. Yes, there are introns and exons in IGH.

2 and 3. You can find sequences and coordinates of IGH genes for human, mouse, pig, teleostei, chondrichthyes, camel, llama, alpaca, caiman, sheep, frog, rat, nonhuman primates, rabbit, platypus, bovine, chicken, dog, Atlantic salmon, rainbow trout, horse, and ring-tailed lemur in the IMGT gene tables:

[IMGT Repertoires][2]

Note: antibodies also contain signal sequences (called L-part on IMGT) that direct the antibody to be displayed on the plasma membrane.

[1]: https://en.wikipedia.org/wiki/V(D)J_recombination
[2]: https://www.imgt.org/IMGTrepertoire/LocusGenes/#F
Question: Hi all,

I have 4 bam files which correspond to 4 runs of one library. I want to merge those 4 files in order to do some variant calling analysis with only one file per sample.

I tried some basic method, with samtools merge, and MergeSamFiles from picard tools, but I am now not sure what they commands do exactly, because at the end I don't find the same number of SNP with one way or the other.

So my question is: What is the method for merge multiple related bam files from the same sample in one? And what is the difference between merge from samtools and MergeSamFiles from picard tools?

I know this is a really basic question and I apologize for that, but I never did this kind of analysis before.

Answer: I am wondering whether the differences you get between using MergeSamFiles and samtools merge may be related to improperly handling read group (RG) information. Each of your BAM files should have a specific RG assigned. This information is taken into account during different steps of variant-calling pipelines: duplicates filtering, variant calling... I would inspect the merged BAM file to see whether RG tags have been properly assigned.
Question: hi, everyone
I am working on calculate the substitution rate of noncoding rnas now. if i have a multiple sequence alignment from three species, the topology of tree is ((1,2),3), how can i use baseml and the REV substitution model to get my desire result? i have tried the example 'brown.nuc' and 'brown.trees' in paml4.7 package, but i cannot make out which number is subtitution rate of that sequences. i have read a lot of papers, but the formula is too difficult to understand, can you give me suggestion?any reply will be appreciated!

Answer: So, it's been a while since I did this type of analysis, but there are two levels of modeling:

1) Modeling substitution rates between two samples (to try and take the probability of multiple substitutions at the same site into consideration). This wikipedia page talks about the models for doing this:

http://en.wikipedia.org/wiki/Models_of_DNA_evolution

2) Using an outgroup (sample 3 in your case) to try and estimate the ancestral state (or, perhaps more precisely, the rate of divergence from an ancestral state). In this particular paper (looking at substitutions at intergenic and synonymous sites), this was done with a relative rate test. See links below for paper and description of strategy:

http://www.plosgenetics.org/article/info%3Adoi%2F10.1371%2Fjournal.pgen.0020163

http://en.wikipedia.org/wiki/Relative_rate_test

From a practical standpoint, I don't remember what specific program outputs look like. However, I know MEGA is typically one of the most popular packages for evolutionary genomics analysis, and I think it has fairly comprehensive documentation:

http://www.megasoftware.net/

If you can get a hold of this textbook, it might also be helpful to you:

http://www.amazon.com/Fundamentals-Molecular-Evolution-Dan-Graur/dp/0878932666
Question: I am trying to download entire dataset for a bioproject using esearch and efetch from the Entrez Utilities.

My syntax is based on syntax posted by @Istvan Albert at https://www.biostars.org/p/111040/#359440, which is

> `esearch -db sra -query PRJNA40075 | efetch --format runinfo | cut
> -d
> ',' -f 1 | grep SRR | head -5 | xargs fastq-dump -X 10 --split-files`

For the BioProject PRJNA269201 I am interested in, slightly truncated syntax as shown below, creates 144 empty files as expected:

esearch -db sra -query PRJNA269201 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | xargs touch

However, when I try the full-length syntax, it behaves differently from what I expected under both scenarios 1 and 2 detailed below:

**Scenario 1**. On head-node of a cluster:

esearch -db sra -query PRJNA269201 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | head -2 | xargs fastq-dump --split-files

one file finished download, but it is 5.5G which is way larger than the 1.2GB I expected based on info at this [link][1] - is this difference because of file compression?! How can I download to a much more compressed version for both storage and downstream RNA-Seq analyses?

-rw-rw-r-- 1 aksrao aksrao 1.1G Jan 19 19:47 SRR1726554_1.fastq

-rw-rw-r-- 1 aksrao aksrao 5.5G Jan 19 19:44 SRR1726553_1.fastq

**Scenario 2**. When I try to submit this as a shell script, the STDERR stream (SLURM queue management on UBUNTU cluster) captures the following error message:

> `2019-01-20T02:28:55 fastq-dump.2.8.2 err: param empty while validating`
> `argument list - expected accession`

This same problem was reported on the original post by user @ bandanaschapagain, but it may not have been answered and resolved, hence I am posting this afresh. Could someone please help me? Thank you!
[1]: https://www.ncbi.nlm.nih.gov/sra/?term=SRR1726553

Answer: I don't think you're doing anything wrong; the first run (SRR1726554) matches what's on SRA (1.1G). I downloaded SRR1726553 myself and also got a .fastq file of 5.6G. It could be that the SRA metadata is wrong; I would contact them and ask for more information (e-mail at ).

You can get a compressed version by calling `--gzip` in your fastq-dump calls. Most aligners will accept gzipped fastq files as input. My full fastq-dump command is:
`fastq-dump $SRA_FILE --outdir $SRA_DIR --gzip --skip-technical --readids --dumpbase --split-files --clip`

I would recommend reading the [Edward's lab fastq-dump article][1] to learn more about some useful options.

For your error message in Scenario 2; I suspect your accession is not getting passed correctly (based on the `expected accession` part.) Maybe add some print calls (like wrapping the fastq-dump part with an `echo` and writing it to a file) in your shell script to see what command it's actually trying execute? The error message seems rather rare, so maybe it's worth asking SRA support about that too.

[1]: https://edwards.sdsu.edu/research/fastq-dump/
Question: I want to normalise my paired-end ChIP-seq data (4 samples treated, 4 samples control) to reads per million (RPM). So far my calculation takes the number of reads mapped to a position, divides it by the total number of mapped reads in the sample, and then multiples this by one million. My question is, should the number of mapped paired-end reads be counted as 1 or 2 in the calculation? I believe Samtools reports the total number of mapped reads as half the actual number of paired-end reads which map?

My current pseudo script for normalising to RPM:

mappedReads = samtools view -c -F 4 input.bam

scalingFactor = 1000000 / mappedReads

bedtools genomecov -ibam input.bam -bg - scale scalingFactor -g chrom.sizes > output.bedGraph

bedGraphToBigWig input.bedGraph chrom.sizes output.bigWig

Answer: You can use [deepTools][1] to do the normalisation that you want. The output is a bigwig file. It is quite fast if you have multiple cores.

bamCoverage -b file.bam --normalizeUsingRPKM -o file.bw

The program first extends the read to match the paired-end length before computing the coverage. I opted to count paired end reads as 2 and not as 1 to avoid a bias when a read is not properly paired which could be a significant fraction.

[1]: http://deeptools.github.io/
Question: I'm working with a [LINCS L1000 dataset][1] that gives the GE of a cell line before and after perturbation by a small molecule. I am using Level 4 data. After loading the .gct file into matlab, I get a matrix of 22268-by-40172 as well as a vector of column_ids and a vector of row_ids.

Using the row ids and the gene metadata txt file included in the download, I know that each row represents a gene.

I can't figure out what a column represents. Obviously, each columns is a single experiment but I can't understand what each id means.

For example, here is a column id "LJP001_BT20_24H_X1_B2_DUO52HI53LO:A03".

So far, I know that "LJP001" refers to LINCS Joint Project and "BT20" refers to the specific cell line. Somewhere, it must contain information about the small molecule used as a pertubagen but I don't know how to interpret this. Any help would be greatly appreciated!

[1]: http://lincsportal.ccs.miami.edu/datasets/#/view/LDS-1203

Answer: To answer my own question.

The column ids for Level 3 and Level 4 data is basically the distil_id. The example I posted

LJP001_BT20_24H_X1_B2_DUO52HI53LO:A03

can be broken into

- the perturbagen group "LJP001"
- the cell line "BT20"
- the brew prefix "LJP001_BT20_24H"
- the plate index "X1_B2_DUO52HI53LO"
- the well index "A03"
- the distil_id "LJP001_BT20_24H_X1_B2_DUO52HI53LO_A03" (note the switch from ':' to '_')

It turns out that the distil_id doesn't contain enough information to identify the perturbagen used. To identify this, you need to use the [LINCS api][1]. Here is [more information][2] about using the LINCS api to query the metadata. I also used [this Coursera video][3] as a reference. Note that the example given in the question doesn't work with the API.

[1]: http://api.lincscloud.org/a2/docs/instinfo
[2]: http://support.lincscloud.org/hc/en-us/articles/202107633-Metadata-Getting-Lists-of-Perturbagens-Cell-Lines-Cell-Types
[3]: https://www.coursera.org/learn/bd2k-lincs/lecture/grLin/introduction-to-lincs-l1000-data
Question: In GSEA, What does it mean enriched gene sets at nominal p-value <1%? and what does it mean also regarding specific phenotype for example (cancer in cancer vs Normal samples)

Answer: Hello

The P-value, or calculated probability, is the probability of finding the observed, or more extreme, results when the null hypothesis of a study question is true. Having a p-value <1% means that is very likely that we will reject the null hypothesis, in your case, that these enriched gene sets are not statistically different between the cancer and the normal samples. In other words, you can say that your enriched gene sets are statistically different between the cancer and the normal samples.

You can watch this video if you have more doubts [StatQuest: P Values, clearly explained][1]

Best

[1]: https://www.youtube.com/watch?v=5Z9OIYA8He8
Question: Hi,

I need to change the affection status in my PLINK files (data1.bed, data1.bim, data1.fam) **to '2' for all samples**. To do this, I was thinking of:

1\. Convert bed to ped using:

plink --bfile data1 --recode --tab --out data2

2\. Copy `data1.fam` as `data2.fam`

3\. Change the 6th column of both the `data2.fam` and `data2.ped` files so that the new affection status is '2' for everything

4\. Convert the modified ped file back to bed using:

plink --ped data2.ped --map data2.map

Do these steps look complete and correct for what I am trying to do?

Thanks

Answer: You could:

1) Change the affected status in the fam file. Smaller file, less steps.

or

2) Use the option `--pheno` to use an alternate phenotype.
Question: Hello Biostars Community,

General question(s) here, specifically in regards to sequencing platforms and some related questions on sequencing, as well?

1. Could using a newer/the newest computational genome/annotation (for example, presently, Ensembl 107 or the newest Gencode version) adversely effect the actual truth of what was sequenced?

2. When sequencing is done through an Illumina machine or other big name company machines, are those sequencing platforms *completely independent* from the genome or DNA/cDNA being sequencing?

3. What happens if, for example, "famous gene ABC" and "low-profile gene XYZ" are found to have different 3' and/or 5' ends by some new discovery, in Illumina, would adapters still link to them to perform those bridge PCR reactions on the flow cell lanes, or would it be that all the data published before on "famous gene ABC" and "low-profile gene XYZ" should be revisited? Or is it like question #2, "completely independent" - are even adapters independent from genes?

I was doing some reading, and I guess the gene sequence really only matters for probe-based sequencing (chips and arrays?). Hopefully this question could be a good resource for others?

Thank you in advance.
- Pratik

Answer: > using a newer/the newest computational genome/annotation (for example,
> presently, Ensembl 107 or the newest Gencode version) adversely effect
> the actual truth of what was sequenced?

No. Annotation represents the current understanding of genome of an organism. Could it have errors? Possibly/more than likely. But those errors get corrected over time (patch releases). Sequence that you got from a run is not going to change. It is independent of annotation/reference you use. Reference used **will** influence your conclusions, so reference will/does play a vital role in final outcome/conclusion.

> are those sequencing platforms completely independent from the genome
> or DNA/cDNA being sequencing?

This will be influenced by limits/characteristics of sequencing technology. e.g. some platforms may not be able to sequence more than a certain number of base homopolymers (most platforms will have some limitations w.r.t this). It may be difficult to get representation of certain areas of genome because they are hard to convert into sequenceable libraries.

> adapters still link to them to perform those bridge PCR reactions on
> the flow cell lanes

Remember that you are adding the adapters to create necessary flowcell compatible ends. The fragments that do not have these ends are not going to bind and will not be sequenced. As long as your fragments have compatible ends (T overlang etc) they will be made into sequence-able libraries.

> gene sequence really only matters for probe-based sequencing

It may matter in case of technologies where you need to be able to unwind the two strands. Sequence may form secondary structures that could be hard to resolve/sequence through.
Question: Hi guys,

I have R programming question: I have more than 1000 genes names `Names`. I want to match the samples containing those genes in `Names` and paste the corresponding gene with the ".AD" extension next to it in the same order to get the `Result` as shown below. Thank you.

```r
Names <- c("cebi", "pithe", "MAPK", "sapiens", "JUNK", "calli", "STR")
samples <- c("MAPK", "JUNK", "STR")
```

`Result`:

"cebi", "pithe", "MAPK", "MAPK.AD", "sapiens", "JUNK", "JUNK.AD", "calli", "STR", "STR.AD"

Answer: ```r
Names <- c("cebi","pithe","MAPK","sapiens","JUNK","calli","STR")
samples <- c("MAPK","JUNK","STR")

Result <- vector()
count <- 1

for(i in 1:length(Names)){
index <- which(samples == Names[i])

if(length(index) == 0){
Result[count] <- Names[i]
count <- count + 1
next
} else{
Result[count] <- Names[i]
count <- count + 1
Result[count] <- paste(Names[i], ".AD", sep="")
count <- count + 1
}
}
```