cannin/biostars_qa huggingface dataset:
Question: Hello,
I am doing a RNA-Seq analysis with R using the limma package. I have RNA-Seq data of different mutant lines of a model organism which have two different phenotypes in comparison with the siblings as control. The samples are like that:

- line 1, phenotype 1, mutant and sibling
- line 2, phenotype 1, mutant and sibling
- line 3, phenotype 2, mutant and sibling
- line 4, phenotype 2, mutant and sibling
- line 5, phenotype 2, mutant and sibling

That means line 1 and 2 have the similar phenotype 1
and line 3, 4, and 5 have the similar phenotype 2
which differ from phenotype 1. Let's say phenotype 1 has more cells than the control and phenotype 2 has less cells than the control.

I did the RNA-Seq analysis for each line, but now I want to compare the lines of one phenotype. I want to find the genes which are different in the lines with the same phenotype in comparison to the siblings. My makeContrast command looks like that:

cont.matrix <- makeContrasts(pheno1 = (line1.mut+line2.mut)-(line1.sib+line2.sib),
pheno2 = (line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib),
levels=design)

Finally I want to find genes (eg. genes involved in the cell cycle) which are in phenotype 1 up and in phenotype 2 down regulated or the other way around.

cont.matrix <- makeContrasts(pheno1vspheno2 = ((line1.mut+line2.mut)-(line1.sib+line2.sib)) -
((line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib)),
levels=design)

I am not sure if a can do a contrast matrix like this. Do you think this will give me the genes I am interested in?

Answer: Let's suppose theres a gene that is 2-fold upregulated in phenotype1 relative to it's siblings and that is also 2-fold upregulated in phenotype2 relative to it's siblings (regardless of the lines).

So this is a gene where the phenotype-induced change in expression is identical.

log2(2) = 1, so we expect line1.mut - line1.sib ~ 1, line2.mut - line2.sib ~ 1, line3.mut - line3.sub ~ 1 and so on.

Now your contrasts:

pheno1 = (line1.mut+line2.mut)-(line1.sib+line2.sib)
= (line1.mut - line1.sib) + (line2.mut - line2.sib)
~ 1 + 1 = 2

pheno2 = (line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib)
= (line3.mut - line3.sib) + (line4.mut - line4.sib) + (line5.mut - line5.sib)
~ 1 + 1 + 1 = 3

and as a result:

pheno1vspheno2 = ((line1.mut+line2.mut)-(line1.sib+line2.sib)) -
((line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib))
= pheno1 - pheno2 ~ 2 - 3 = -1

So, although this gene was induced to the exact-same-extent in mutant vs sib in each of your individual lines, there appears to be a two-fold-lower induction in phenotype1 relative to phenotype2. So something is wrong with the maths underlying your contrasts.

To prevent this, you should normalise your contrasts by the number of lines:

pheno1 = (line1.mut+line2.mut) / 2 - (line1.sib+line2.sib) / 2

pheno2 = (line3.mut+line4.mut+line5.mut) / 3 - (line3.sib+line4.sib+line5.sib) / 3

pheno1vspheno2 = ((line1.mut+line2.mut)-(line1.sib+line2.sib)) / 2 -
((line3.mut+line4.mut+line5.mut)-(line3.sib+line4.sib+line5.sib)) / 3
Question: Hi everyone,

I am using the TCGA portal to get mRNA expression data for various cancer studies (e.g. lung, liver, thyroid etc). I have two questions about the data:

1. Some cancer studies on TCGA have "mRNA expression RNASeq V2 RSEM" values & corresponding "z-scores". I am confused as to what the "mRNA expression z-Scores (RNA Seq V2 RSEM)" data constitutes of. How are the z-scores calculated and what do they represent?
2. We have been on a lookout for control dataset for the cancer studies on TCGA. Does anyone know of a good place where you can find control dataset for tissues like Lung, Liver, Thyroid etc. (basically all the fore-gut tissues)? We are working with control data from [GTEx][1] but they have RPKM values & TCGA has RSEM/RSEM z-scored values, so we have to do a lot of scaling/normalization/transformation to compare these disparate datasets. We would like to know if there is any mRNA expression data (obtained via RNASeq V2 RSEM) for controls.

Thanks!

UPDATE: I have posted the second part as a separate question [here][2].

[1]: http://www.gtexportal.org/home/
[2]: https://www.biostars.org/p/107594/

Answer: I will attempt to answer part 1 of your question:

A z-score for a sample indicates the number of standard deviations away from the mean of expression in the reference. The formula is:

z = (expression in tumor sample - mean expression in reference sample) / standard deviation of expression in reference sample

TCGA states:

> For mRNA and microRNA expression data, we typically compute the relative expression of an individual gene and tumor to the gene's expression distribution in a reference population. **That reference population is either all tumors that are diploid for the gene in question, or, when available, normal adjacent tissue.**

(It is not always clear what the cell of origin of a tumor is, so the mRNA expression in normal adjacent tissue can sometimes be misleading, which is why expression is sometimes compared within the set of tumors only).

As for part 2, CPM (count per million) data for each gene and sample would be ideal for the cross-sample comparisons, but I am not sure where you could get such data. Maybe you should post this as a separate question?
Question: How much does it cost to align a 8 lanes of a HiSeq 2000 flowcell paired-end run consisting of 6 billion 100bp human genomic reads (600Gbp output) using the Amazon EC2 and associated storage?

Assume using an aligner such as BWA with default parameters, aligned against a human reference. Please compute the cost and time from upload of FASTQ files to download of BAM files inclusive.

"It depends" is not a useful answer - if there is a factor such as instance-type or EBS vs S3 please select a sensible option and provide a quote.

Answer: Konrad Karczewski of stormseq said in a recent tweet that "a full genome (30X coverage) is about $30 or so, exome (80X coverage) is around $2-3." I think this is from fastq->vcf using BWA and GATK.

Assuming an 80X exome is about 40M reads, a $2.50 exome comes out to about $6.25/hmmvr (hundred million mapped variant reads - my awesome new metric), or about $375 for a HiSeq 2k flowcell. Granted this is going all the way to VCF but that number is in line with what I've heard elsewhere.

I think the BWA->GATK pipeline is commoditized enough that providers should be using these $/hmmvr on scienceexchange. Then labs will have something to compare it with, if they don't want to try AWS themselves.
Question: Hi.

I made vcf files using UnifiedGenotyper in GATK and found that many of mutants had 1/2 at genotype (GT) field.

As I understand, 0/0 = homo ref, 0/1 = hetero, 1/1 = homo alt.

Does anyone know what exactly 1/2 mean?

Thank you!

Answer: There are apparently two alternate alleles at this position. So if the reference is `C`, then perhaps you have an `A` and a `G`.
Question: I have to calculate the Shannon entropy for a given list of sequences in a fasta file. Recently I came across a program called rnaentropy which solves my problem but the issue is I cannot run sequences in batch. In the documentation the input is given as a sequence of strings. The parameters in the documentation is as follows:

./RNAentropy "sequence" -s sequence -t temperature -e energyModel(default 2004) -d delta_Temp (compute structural entropy using <E> = RT² * d/dT ln(Z(T)) ) -c (centered) -z energy_is_zero [0|1] (default 0) -v (verbose, extended output)
./RNAentropy -h (detailed help)

Input paramters are:
<sequence> : If FIRST argument is a valid nucleotide sequence it will be used input
-s <sequence> : Alternative flag for input sequence
-t <temperature> : Temperature in ºC (default is 17ºC)
-e <energy_model> : Thermodynamic energy model used. Valid values for energy model are 1999, 2004 and 2007 (resp. Turner'99, Turner'04 and Andronescu '07) (default is 2004)
-d <delta_T> : Temperature variation used for estimating d/dT ln(Z(T))
NOTE: If this parameter is provided, H is computed estimating <E> = RT² * d/dT ln(Z(T))
-c : (Use only in combination of -d) Use the centered version for estimating <E> = RT² * d/dT ln(Z(T))
-v : Output includes method for computing H and the name of the ouput parameters
-z <energy_is_zero>: [0|1] If value is 1, energies are set to 0, ouput is structural entropy for the uniform case

Can anyone please suggest me a way to do it for a multi-fasta? Is there a way to iterate in a shell script?

Answer: I wouuld use a Perl (or Python script) to parse the multifasta and run your tool, like:

#!/usr/bin/perl

use strict;
use warnings;

$/="\n>"; # read fasta sequences in blocks
while (<>) {
s/>//g;
my ($id, @seq) = split (/\n/, $_);
my $seq = join "", @seq; # rebuild the sequence, one single line
print "$id\n"; # to indicate which sequence is being analyzed
system("./RNAentropy -s $seq"); # run the command, add params if needed
}

then you can run as `perl analyzeSeq.pl < multifasta_in > output_file`
Question: Are there any programs available (Python, Perl, whatever), that can take a degenerate nucleotide sequence and translate it into its multiple possible oligos?

Any help would be appreciated.

Answer: This program should do the job:

```
#include
#include
#include
#define RUN(base) copy[index]=base; recurs(seq,copy,index+1,len)

static void recurs(const char *seq,char* copy,int index,const int len)
{
if(index==len)
{
fwrite(copy,sizeof(char),len,stdout);
fputc('\n',stdout);
}
else
{
switch(toupper(seq[index]))
{
case 'A':case 'T': case 'G': case 'C': RUN(seq[index]); break;
case 'N':RUN('A');RUN('T');RUN('G');RUN('C');break;
case 'W':RUN('A');RUN('T');break;
case 'S':RUN('G');RUN('C');break;
case 'B':RUN('T');RUN('G');RUN('C');break;
case 'D':RUN('A');RUN('T');RUN('G');break;
case 'H':RUN('A');RUN('T');RUN('C');break;
case 'V':RUN('A');RUN('G');RUN('C');break;
case 'K':RUN('G');RUN('T');break;
case 'M':RUN('A');RUN('C');break;
case 'R':RUN('A');RUN('G');break;
case 'Y':RUN('C');RUN('T');break;
default: fprintf(stderr,"Bad base in %s (%c)\n",seq,seq[index]); exit(EXIT_FAILURE);break;
}
}
}

int main(int argc,char** argv)
{
char* seq;
int len,i;
if(argc!=2)
{
fprintf(stderr,"Usage : %s ",argv[0]);
return EXIT_FAILURE;
}
seq=argv[1];
len=strlen(seq);
char* copy=malloc((len+1)*sizeof(char));
if(copy==NULL)
{
fprintf(stderr,"Out of memory\n");
exit(EXIT_FAILURE);
}
copy[len]='\0';
recurs(seq,copy,0,len);
free(copy);
return 0;
}
```

**Compilation**:

gcc -O3 -o prg prg.c

**test**

```
> ./prg ATGCTGATCGAGCTANATCGATCGGACTACY
ATGCTGATCGAGCTAAATCGATCGGACTACC
ATGCTGATCGAGCTAAATCGATCGGACTACT
ATGCTGATCGAGCTATATCGATCGGACTACC
ATGCTGATCGAGCTATATCGATCGGACTACT
ATGCTGATCGAGCTAGATCGATCGGACTACC
ATGCTGATCGAGCTAGATCGATCGGACTACT
ATGCTGATCGAGCTACATCGATCGGACTACC
ATGCTGATCGAGCTACATCGATCGGACTACT
```
Question: Dear all,

I would like to align some reads to a reference using STAR. The following command works perfectly:

STAR --genomeDir output/spiking/index/star --readFilesIn reads.fastq --outFileNamePrefix outputFolder --runThreadN 8 > message.txt

However, the following command does not work:

STAR --genomeDir output/spiking/index/star --readFilesIn reads.fastq.gz --outFileNamePrefix outputFolder --runThreadN 8 > message.txt

Error: the read ID should start with @ or >

Presumably, STAR expects a FASTQ file and not a fastq.gz file. Does anybody know how to get round this in an efficient way?

Thanks.

C.

Answer: add `--readFilesCommand zcat` in your commandline
Question: Do you know any public scientific SQL server ?

for example, I would cite:

UCSC http://genome.ucsc.edu/FAQ/FAQdownloads#download29

ENSEMBL http://uswest.ensembl.org/info/data/mysql.html

GO http://www.geneontology.org/GO.database.shtml#mirrors

(I'll give a +1 to each correct answer)

Answer: PublicHouse - uses the BioWarehouse system; requires user registration.
Question: Hi, I want to use bedops to analyze some of the .bed files from ChIP-seq, but one of the files I can't go through the analysis, even can't sort-bed. It keeps saying the BED row length exceeds capacity at line 1. It doesn't help even after I deleted the first row of my .bed file.

```
unknownb8f6b1106ced:chip-seq fuxiaoyong$ sort-bed era_inpegf_mcf7_n3_hg18_f20_nr.bed
BED row length exceeds capacity at line 1 in era_inpegf_mcf7_n3_hg18_f20_nr.bed.
Check that you have unix newlines (cat -A) or increase TOKENS_MAX_LENGTH in BEDOPS.Constants.hpp and recompile BEDOPS.
```

I am very new for using tools to analyze NGS data, even for the Mac OS X system. So, please help me and your detailed explanation will be greatly appreciated!!

Answer: Hello,

As discussed, the problem was owing to Excel messing up your BED file. To delete blank rows from your BED file, use the following command:

sed '/^$/d' era_inpegf_mcf7_n3_hg18_f20_nr_NoBlanks.bed

Or, you could also manually remove blank lines (a tedious process) by opening the file in TextWrangler.