cannin/biostars_qa huggingface dataset:
Question: Hello everyone,

I'm working with data that was performed using Agilent-014850 Whole Human Genome Microarray 4x44K G4112F arrays and the One Color Quick Amp Labeling Kit (Agilent Technologies, Santa Clara, USA).

I'm trying to set up a code in R:

#Load limma
library(limma)

#Set-up
targetinfo <- readTargets("Targets.txt",row.names="FileName",sep="")

#Read files
project <- read.maimages(targetinfo,source="agilent", green.only=TRUE)

#Background correction
project.bgc <- backgroundCorrect(project, method="normexp", offset=16)

#Normalize the data with the 'quantile' method for 1-color
project.NormData <-normalizeBetweenArrays(project.bgc,method="quantile")

#Create the study design and comparison model
design <- paste(targetinfo$Target, sep="")
design <- factor(design)
comparisonmodel <- model.matrix(~0+design)
colnames(comparisonmodel) <- levels(design)
#Checking the experimental design
design
comparisonmodel

project.fit <- lmFit(project.NormData, comparisonmodel)

> #When I used > project.fit <- lmFit(project.NormData$M, comparisonmodel)
> Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x), :
> 'data' must be of a vector type, was 'NULL'

> #So, I used without a $M
> project.fit <- lmFit(project.NormData,comparisonmodel)

#Applying the empirical Bayes method
project.fit.eBayes <- eBayes(project.fit)
names(project.fit.eBayes)

#Make individual contrasts and fit the new model
CaseControl <- makeContrasts(CaseControl="D0A-Control", levels=comparisonmodel)
CaseControl.fitmodel <- contrasts.fit(project.fit.eBayes, CaseControl)
CaseControl.fitmodel.eBayes <- eBayes(CaseControl.fitmodel)

#Filtering Results
nrow(topTable(CaseControl.fitmodel.eBayes, coef="CaseControl", number=99999, lfc=2))
probeset.list <- topTable(CaseControl.fitmodel.eBayes, "CaseControl", number=99999, adjust.method="BH", sort.by="P", lfc=2)

#To save results
write.table(probeset.list, "results.txt", sep="\t", quote=FALSE)

Well, now it seems correct, and you guys what do you think?

Best regards,

Leite

Answer: You have to pass a comparison model to `lmFit()`. Also, for Agilent processing, the expression values are eventually stored in the 'M' or 'E' variables of your object, and not accessed via the `exprs()` function. The choice of M or E depends on whether it is 2- or 1-channel / colour array.

comparisonmodel <- model.matrix(~0+design)
colnames(comparisonmodel) <- levels(design)
comparisonmodel

fit <- lmFit(project.NormData$M, comparisonmodel)
Question: Hi there !

I am working on Ribosome Profiling data, and I want to use Python3 multithreading to speed up my program.

To describe the scope of my question:

I use python scripts that take as input alignment files like:

python3 python_script.py <(samtools view file.bam)

and I take it in charge in the script with the fileinput module (**import fileinput**). It works like a charm !

So related to this, do you think that it is safe to call 3 python scripts at the same time using the same file.bam (using the samtools view command) ? For instance:

python3 script1.py <(samtools view file.bam) &
python3 script2.py <(samtools view file.bam) &
python3 script3.py <(samtools view file.bam) &

I am not sure: if the samtools view command is called on the file.bam, can I safely read again this file with an other samtools view command ? The output is not a problem, every scripts will write a different file. I gave the example with 3 scripts, but I will certainly use 3 functions.

Any kind of help of suggestion will be highly appreciated.

Wocka

Answer: You can call samtools (or any other program) any number of times on the same BAM file without issue. You can also use pysam inside of python to open the same files multiple times within the same script without issue.
Question: This question is similar to this one (Get rs number based on position).

I have a text file with SNPs in the chr:position format

10:71086
10:72876
10:75794

I was wondering if there is an R package (or perhaps one in Bioconductor) that can take these as input values and provide the SNP rs number associate with them????

Answer: I think JC has already pointed out a nice resource. The other option would be

Get rs number based on position
Question: Hi All,

I have a table with six columns (in the following).

![enter image description here][1]

[1]: http://uupload.ir/files/xm5g_taible.png

So, I want to delete SNPs that are less than number 10 in the SNP column.

What is the best idea?

Best Regard

Mostafa

Answer: Read the file, and filter out any row where SNP column is more than or equal to 10:

data <- read.table("SNP.txt", header = TRUE)
data <- data[ data$SNP >= 10, ])

To exclude 1st BP column, and to add 20KB to BP column:

data <- data[ , -3 ]
data$BP <- data$BP + 20000
Question: Hi everyone!

I am trying to find the best way to make 2 boxplot for a specific gene from data found in a row for a subset of columns within data frame x.

x dimensions are 634 by 128 columns

Each row is specific to a gene,

Column 1 has gene name, and I want to say look at gene in row#1

columns 2:48 data I want to include in one boxplot

columns 49:128 data I want to include in another boxplot

data frame looks something like this

```
gene accepted_hits_x1.bam accepted_hits_x1.bam etc....
1 AARS1 -6 0 etc....
```

I also want to be able to see each data point that makes up the boxplot plotted in the plot

I am having a problem:

I am running into the problem where my data (residual from mean ... meaning x value - mean) is a series of positive and negative values and it appears that with this plot it is excluding these negative values...

```r
data <- unlist(subset(datavr, gene =="IGF1R", select=2:128))

news <- data.frame(data=data, factor=c(rep(1,47), rep(2,80)))
news$data <- (log10(as.numeric(news$data)) + 1)

g <- ggplot(data=news, aes(x=as.factor(factor), y=data))

g + geom_boxplot() + geom_point(color="purple", size=3) + xlab("A38-41 A38-5 ") + ylab("log10(Residual from Mean)+1") + ggtitle("IGF1R inside region") + theme(plot.title = element_text(face="bold"))
```

The problem is that it keeps giving me error saying that:

```r
Removed 110 rows containing missing values (geom_point)
```

This could be that these values are negative so taking the log10(value)+1?

Answer: To make a boxplot with base graphics in R, you need to create a "factor" vector, which indicates which category each of your data points belong to:

```r
f <- factor(c(rep("Group 1", 47), rep("Group 2", 80)))
```

Then call "boxplot" with the factor and data as arguments:

```r
boxplot(f, as.numeric(dat[1, 2:128]))
```

Edit: Actually creating the factor is not even necessary in this case. You can just list the two data vectors as multiple arguments:

```r
boxplot(as.numeric(dat[1,2:48]), as.numeric(dat[1,49:128]))
```
Question: Hi,
I recently found very helpful http://cancergenome.broadinstitute.org/ portal to find out which gene is mutated in which cancer, frequencies and vice versa. But I can't find out how to determine frequencies for the combinations of genes. For example, based on data from this portal PTEN is mutated in 7% of all cancers and in 30% of gliobastomas. Rb is mutated in 3% of all patients and in 8% of glioblastomas. But what will be the easiest way for non bioinformatician to find out the frequency of patients with both Pten and Rb mutation in general population, in glioblastomas and in other cancers?

Thanks for your help!
Andrey

Answer: You can get at this kind of information using the cBioPortal for Cancer Genomics.

To get a relevant visualization:

Go to http://www.cbioportal.org/public-portal/

Enter your gene names using their official HUGO designations: PTEN, RB1. One per line.

If you are only interested in mutations, select the 'Only Mutation' option

Hit the submit button

Hover your mouse over the bar for your cancer type of interest and select 'click to view details'

This will bring you to an 'OncoPrint' that shows which cases have mutations in: PTEN, RB1, both, neither

Also check out the additional tabs of this view to explore the data for these two genes further. You can also refine your query using the 'Modify Query' section.

Here is an example query for: RB1 and PTEN in GBM

In this example, 31% have PTEN mutations, 9% have RB1 mutations, and 5% have both. According to the 'mutual exclusivity' tab there is a highly significant tendency tendency toward co-occurrence of mutations in these mutations (odds ratio: 3.875, p-value: 0.001463)
Question: Hello,

My question is how to predict one connected protein structures from two separate known structures?

It should be different from the protein-protein docking method, which is non-covalent interactions. My question is that the two proteins should "contently" link together.

Is that possible first manually link them (two 3D structure) together using some software, and then do the energy minimization?

Any suggestions are appreciated!

Answer: What you want to do is possible, but can't be done reliably without external restraints. Protein modeling typically means satisfying a number of angle and distance restraints. If you add a linker connecting two domains of known structures, there will be no restraints specifying their relative orientation unless you add them manually. Energy minimization will help you only if you are already close to the optimal domain interaction, but it isn't meant for making proteins move around a lot.

There are many recent methods predicting intra-protein contacts between protein residues based on co-evolution. Those approaches can be extended to inter-protein contacts as well, and may provide some restraints for your modeling procedure.

- https://academic.oup.com/nar/article/46/W1/W432/5001161
- https://link.springer.com/protocol/10.1007%2F978-1-4939-9873-9_6
Question: Hello,

My name is Maria. I'm a master's student from Estonia. I want to ask some basic (okay, maybe a little bit silly questions). I'm starting to work with 1000Genomes data and I've never worked with genome sequences before.

I want to download sub-sequences of a genome. The instruction says to indicate it like 1:1-50000. I understand that 1 in front of : refers to chromosome number, is that correct? And 1-50000 would be the first 50000 nucleotides?

Are the genomes of different people of different lengths due to copy number variations? Does the sequencing according to a reference genome take account of these differences or would all the genomes in 1000Genomes be of the same length as they are aligned to a reference genome?

What if I wanted to obtain specific genes from the sequences? Is there any tool to do that?

Thank you!

Answer: > I want to download sub-sequences of a genome. The instruction says to indicate it like 1:1-50000. I understand that 1 in front of : refers to chromosome number, is that correct? And 1-50000 would be the first 50000 nucleotides?

Yes, you are correct

> Are the genomes of different people of different lengths due to copy number variations?

Yes. Indels will also contribute to differences in genome lengths.

> Does the sequencing according to a reference genome take account of these differences or would all the genomes in 1000Genomes be of the same length as they are aligned to a reference genome?

Reads from different genomic samples are aligned to the same reference genome so that multiple genomes can be easily compared to each other for the presence/absence of a genomic variant. Otherwise it would be tough to carry out any comparisons. In short, you can say that all the genomes in 1000 Genomes are of same length w.r.t to the coordinate location of a variant/gene. Though bam files have enough information to predict copy number variants and identify insertions and deletions differing between individuals.

> What if I wanted to obtain specific genes from the sequences? Is there any tool to do that?

You can download coordinates of your gene of interest from Ensembl (gtf file) or UCSC genome browser (gtf/bed) and then use those coordinates (for e.g. chr2:100000-1020000) to fetch the reads overlapping that region from the bam file. You can use samtools view function.
Question: I would like to add read group info (`-R`) during the mapping/alignment stage as part of my variant calling gatk pipeline.

I am doing something like this:

bwa mem \
-M \
-t 8 \
-v 3 \
-R <(sh a-illumina-read-group.sh $1) \
"$path_dr_bwaindex_genome" \
$1 $2

where the read group string is generated automatically from the fastq file by the shell script `a-illumina-read-group.sh`. It produces a string like:

`'@RG\tID:ST-E00215_230_HJ3FMALXX_2\tSM:ST-E00215_230_HJ3FMALXX_2_ATCACG\tLB:ATCACG\tPL:ILLUMINA'`

but, when I run bwa, it fails with this error:
`[E::bwa_set_rg] the read group line is not started with @RG`

I have tried excluding the single quotes (`''`) around read group info, but that didn't change anything. I also tried variations in how the variable is passed.

`-R=$(sh a-illumina-read-group.sh $1) \`

`-R=$(echo $(sh a-illumina-read-group.sh $1)) \`

Just as a test, I tried:

`rg='@RG\tID:ST-E00215_230_HJ3FMALXX_2\tSM:ST-E00215_230_HJ3FMALXX_2_ATCACG\tLB:ATCACG\tPL:ILLUMINA'`

`-R <(echo $rg) \`

But, they all produce the same error. I would appreciate any solutions to this issue.

I understand that I can add read groups later on using PicardTools `AddOrReplaceReadGroup`, but I thought it might be convenient doing it here in one step.

Thanks.

Answer: What worked for me was to read the read group information from the fastq file during the mapping run. My read name in the fastq file looks like this: `@ST-E00274:188:H3JWNCCXY:4:1101:5142:1221 1:N:0:NTTGTA`.

The bash file looks like this:

#!/bin/bash

header=$(zcat $1 | head -n 1)
id=$(echo $header | head -n 1 | cut -f 1-4 -d":" | sed 's/@//' | sed 's/:/_/g')
sm=$(echo $header | head -n 1 | grep -Eo "[ATGCN]+$")
echo "Read Group @RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA"

bwa mem \
-M \
-t 8 \
-v 3 \
-R $(echo "@RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA") \
"$path_bwaindex_genome" \
$1 $2 | samblaster -M | samtools fixmate - - | samtools sort -O bam -o "mapped-bwa.bam"

And the bash file is run as

bwa-mapper.sh read_1.fq.gz read_2.fq.gz

You can remove this part (`| samblaster -M | samtools fixmate - - | samtools sort -O bam -o`) and replace with `>` if you don't need it. `samblaster` marks duplicates like Picard. I don't remember what `fixmate` does. The sort part sorts the SAM file and generates an output BAM.
Question: Hi, I'm new to both proteomics and biostar and would like to start with maybe a simple, but a more general question.

I know that for RNA-Seq, people just assume normality to the data mostly, because to achieve normality with such a small number of samples is always risky and difficult. I was wondering how this is handled in the proteomics field. What kind of tests one might use to analyze the data, when a t-test is not applicable due to the non-normal behavior of the data?

What alternatives do I have there?

I saw that some people are using LIMMA. Can LIMMA handle non-normal distributed data?

thanks

G.

Answer: When using `limma` in proteomics/metabolomics, the data are usually log-transformed first, which is usually makes it normal/more homoscedastic (or assumed to be so). (Example [paper][1]).

By the way, people do not usually assume normality for RNA-seq data. They usually model raw counts with negative binomials (e.g. `DESeq2`, `edgeR`) or do some other transformations (e.g. `limma-voom`).

[1]: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5418186/
Question: Hi there,
I'm new in bioinformatics tools and I need help.
I tell you what I have done :)

I have an assembly from Trinity and I wanted to know the CG content of each assembly so I generated a file with this information and I imported to R.

Then, I calculated the CG content of all my contigs in R by doing:

## load data
data= read.table("CG_content_contig.txt", header = T)

## Create a new variable called CG
data$CG <- (data$C+data$G)/(data$A+data$C+data$G+data$T)

## Create a file with CG content > 23%
CG_content_more_than_23 <- data[which(data$CG > '0.23'), ]

## Create a file only with the names of the contigs
name_contigs <- CG_content_more_than_23$chr

## Download .txt file.
write.table(name_contigs, file="name_contigs.txt", row.names=FALSE, sep='\t')

As you see, I have created a file with those contigs with an CG content higher than 23%. I have uploaded this file (name_contigs.txt) to my server in order to work with it. It looks like this:

TRINITY_DN134693_c0_g1_i1
TRINITY_DN109669_c0_g1_i1
TRINITY_DN109679_c0_g1_i1
TRINITY_DN114999_c0_g1_i1
TRINITY_DN114910_c0_g1_i1

I have a .fq file that looks like this:

>TRINITY_DN134617_c0_g1_i1
AATAAAAATAAATAAAAATCAATAAAAATATTATAATACAATATAATATAAAATAATATAAAAATTCTACAATAAGAATAAAGTATAATTTTTTAGATTATAAGAGGATATGTTAATACATAGTATTCTGTTTGTTATTGTAGAAAAAACATACAGAAACTTTTTGTATATATAGTCTCATTTTATATATATAAATAAAAATGAACATTAATGAAATGAAATTAAGAGTCGTTTTATTAAAAATAGCTATAAAAAATAACAACA

>TRINITY_DN134643_c0_g1_i1
GCATGGTAGTAAAGTATAATGACATAGCAAAAATATTTAAAATAAAAAAAAATTACTATTATAATTTTTTCTGTATAACATAAACGTTTTTAATGATATTATATTAATTACATATAAAAATAGCATAATAAAAATATTTAGTTATAAAATTTATTATTTTATTTTTTTTTTTTTGTTATATACTTTCTCAGAACATTAATTTGTCATCAGTTCTATTATATTGATAAACTATTCAATTGCTTTAATA

What I want to do is to keep only the contigs that are NOT in the .txt file.

Maybe I should create a pipeline with grep command?

Answer: You should use `faSomeRecords` [utility from Jim Kent][1] (UCSC, linux binary linked). Add execute permissions after download `chmod u+x faSomeRecords`.

faSomeRecords - Extract multiple fa records
usage:
faSomeRecords in.fa listFile out.fa
options:
-exclude - output sequences not in the list file.

[1]: http://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faSomeRecords
Question: Hi everyone!

Sorry for this noob question, but I just wanted to ask, what is the best practice in the field for publishing RNA sequencing data analysis scripts? Does it usually happen before the publication, or after? Currently I kept the repo as private and the paper is not published yet.

Answer: I would strongly recommend publishing them before the paper. That way, you can get your code peer reviewed before you submit your article. That may improve your code. You may have to arrange reviewers yourself for this (maybe an idea for a new BioStars section?).
Question: Hello, I'm using deeptools computeMatrix/plotHeatmap, and am getting a large amount of missing data (i.e. black in the graph).

Commands used
----------
The call to computeMatrix takes in ChIPseq peaks, and corresponding signal in bigWig format:

**Command 1: computeMatrix**

computeMatrix reference-point --referencePoint center -S file1.bw file2.bw file3.bw file4.bw -R peaks.bed -a 500 -b 500 -o matrix.gz

*Note: actual filenames, -p option are omitted*

*Note: --binSize is omitted or 1*

This is followed by:

**Command 2: plotHeatmap**

plotHeatmap -m $infile -o $ofile --outFileSortedRegions ${ofile/.$outformat/.bed} --sortUsing mean --outFileNameMatrix ${ofile/.$outformat/.mat} --zMin 0

*Where: infile=matrix.gz (output from first command); ofile=matrix.$outformat; outformat typically equals "png"*

Output
----------
I have verified that the data are indeed not missing, in at least the first few examples. See the highlighted area of the attached graph (red box).

![plotHeatmap screenshot][1]

Here is a chunk of the corresponding data stored via the --outFileNameMatrix option:

> 1.256 1.256 1.256 1.256 1.256 1.256 1.256 1.256 1.256 1.256 1.346 1.346 1.346

I've validated these numbers by looking specifically at this region of the input file (equivalent to ***file1.bw*** in **Command 1**)

#converted to wiggle
track type=wiggle_0
fixedStep chrom=chrII start=2378124 span=1 step=1
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.25641036
1.34615386
1.34615386
1.34615386

*Note: this matches the above matrix output when --binSize is set to 1 for **Command 1: computeMatrix**; the default setting (10) gives the expected results.*

This value range should be printed as yellow according to the figure legend (not in screenshot), and represents the middle range of the data being plotted.

Also, adding the computeMatrix parameter --missingDataAsZero, you would expect all of that black to become red, right?
That is not reflected in this plot:

![plotHeatmap with black missing data versus zeroes][2]

## Question ##

Does anyone know what might be going on here? Am I running the commands incorrectly (newb)?

[1]: https://s8.postimg.cc/xjbyknr3p/plot_Heatmap_screenshot.png
[2]: https://s8.postimg.cc/kqe8xoyxx/Na_Ns_As_Black_Versus_Na_Ns_As_Zero.png

Answer: Every time this issue is brought up it turns out to be due to the interpolation done in order to make a large matrix fit into a much smaller number of pixels. Try doing one of the following:

1. Make the image larger so there's less interpolation.
2. Increasing the DPI
3. Changing the `--interpolationMethod` option

I have yet to see an example where that doesn't produce the results you're expecting. At the end of the day this comes down to how best to compress a lot of information into a small image. At some point you have to throw information away since there aren't enough pixels to represent everything, so no particular method will always produce the most desirable results.
Question: hi,

I have a long list of different miRNAs like below

hsa-miR-7641
bta-miR-2904
hhi-miR-7641
hsa-miR-4454
bta-miR-2478
efu-let-7c
mmu-miR-6240
hsa-miR-7704

is there any way to avoid searching one by one in mirbase and copy pasting their mature sequences and taking their sequences at the same time?

Answer: It looks like you can download sequence files:
http://www.mirbase.org/ftp.shtml

"Fasta format sequences of all miRNA hairpins" & "Fasta format sequences of all mature miRNA sequences"

When you have those files, a simple `grep -f ` would do the trick
Question: I'm using the equation

Num_Reads * Avg_Read_Length / Genome Size

To calculate coverage. Here are my questions

1. Should I only consider mapped bases in the Query to calculate the read length? Specifically, `M`,`=` tags in the CIGAR string? For example, if a read has `60M10S` for a CIGAR string, would the read length be 60 since 10 were not aligned properly?

2. Given 1. why is it that my calculation of coverage differs from `samtools mpileup`? Wouldn't the average of column 4 in mpileup output equal the average coverage?

Thanks for any help.

Answer: 1) yes M , = and even X

2) yes , ignore the clip

3) did you consider some options like:

-Q, --min-BQ INT skip bases with baseQ/BAQ smaller than INT [13]
--ff, --excl-flags STR|INT filter flags: skip reads with mask bits set
[UNMAP,SECONDARY,QCFAIL,DUP]

what's the output of `samtools depth`
Question: Hi all,

I think GATK is a great toolbox. There a quite a few steps involved and I was wondering on the impact and importance of joint genotyping - in particular when working with very small sample sizes (around 10 -15 samples).
I read that it can lead to a loss on unique snps, which we would be particularly interested in. Has anymore experience with this and can tell me about the effects on the outcome? In particular in quantitative terms, do I gain many more overlapping SNPs - do I loose a lot of unique SNPs? Or is the impact negligible?

Thanks for your input!

Answer: There is a great answer provided previously: https://www.biostars.org/p/10926/

To summarize:

> joint calling has biased false negative rate: it does better if a SNP
> is shared between samples but worse if it is a singleton