cannin/biostars_qa huggingface dataset:
Question: hi all

I am working with somatic mutations in prostate, I downloaded the MAF files + clinical data from TCGA and I installed MuSiC to work with them...however, I need bam or wig files. since I have no access to the bam files, I would like to get hold of the Wiggle files but I am not really sure where to find them/which ones to use...

also, wig files are still pretty big guys (i.e. wiggle files for some ChIPSeq data I worked with are >200Mb each...) ...so if I need to download 200Mb x 400 patients that's gonna be pretty challenging..any suggestions how to circumvent that?

any help would be great.

thanks

Answer: Wiggle files contain coverage status per bp, that some TCGA GSCs like Broad upload within their MAF archives. So for TCGA PRAD, you'll find them [here][1].

~2 years ago, I had tried to systematically collect coverage data across TCGA exome-seq. It doesn't include a lot of the newer data. But all my notes are now belong to you: http://dx.doi.org/10.7303/syn1695394

[1]: https://tcga-data.nci.nih.gov/tcgafiles/ftp_auth/distro_ftpusers/anonymous/tumor/prad/gsc/broad.mit.edu/illuminaga_dnaseq_curated/mutations/broad.mit.edu_PRAD.IlluminaGA_DNASeq_curated.Level_2.1.4.0/
Question: When using the tab output format, it causes multiple outputs of the same hit to show, but with different %identities and other attributes.

When using the standard output, the results are what I find on the web based version of the blastn tool.

How can I create a tabular form of results without getting all these duplicate hits?

Answer: Those are not duplicates, those are partial alignments.
Question: how do I convert fa files to bed format?

Answer: This is a good, though slightly misguided question. If you want to make a BED file from a FASTA sequence, you might do something like this:

Find your FASTA sequence. We'll use human hemoglobin beta as an example.

Use BLAST to align your sequence to the human reference genome.

In the alignment example above, you would pick the genomic alignment, not transcript, and choose "subject" start and end positions. Note that this example has 3 exons, so your BED file start would be 5186957 and end would be 5188159.

Create a BED file using the chromosome (11 in above example), chromstart (5186957), chromend (5188159), and a name for your gene (HBB). You can then upload this as a custom track in the UCSC genome browser.

You will probably want to add introns to your gene structure, so at this point you can the GenePred format instead of BED. This allows you to specify the strand, as well as introns, exons, and transcription factor binding sites.
Question: Hello,

I am new to bioinformatics. I am learning about SNPs and trying to use vcftool. In vcftool if i use the --012 option it generates a .012 file with 0,1,2 in it. Where "*Genotypes are represented as 0, 1 and 2, where the number represent that number of non-reference alleles.*"(Ref: Vcftool manual). So If for a snp position the reference allele is C the possible combinations are CC,CT,TC,TT where CC=0, CT/TC=1, TT=2 are the encodings.

My first question is what does this CC mean ?? Does this mean at the particular position there is neocleotide C and it is paring with another neocleotide C on the other strand?

If that so it is not following the [base-pair rule][1] . Is it ?

Finally the generated files have very large number of 0's and 2's than that of number of 1's. Should not CT/TC be the common combination that CC and TT.

If I am understanding completely wrong some helpful links will be appreciated.

[1]: https://en.wikipedia.org/wiki/Chargaff%27s_rules

Thanks In Advance.

Answer: CC means homozygous for C. CT/TC would mean heterozygous and so on.
Question: I have some fastq files that contain a sequence length distribution that are not uniform. Most sequences are trimmed at 75bp with ~3/8ths at 74bp ~1/4th at 73bp and trailing off all the way to the 34bp length. This closely resembles an exponential curve. Anyway, I am curious if the varying lengths will cause problems when attempting to use TopHat2 to align to a reference genome.

Answer: Tophat2 supports variable read-length. So using Tophat2 won't be an issue. Personally, I don't use reads with length less than 40 nt. I don't think reads with variable read lengths will cause any problem with the interpretation of results generated from Topaht2.
Question: I don't understand what differential expression is, and an article I read mentioned it. I can't find anything to actually answer my question. Can anyone help me better understand what this means?

The link to the article is: https://tcr.amegroups.com/article/view/957

Answer: Very basic summary:

Let's say that you have two samples from the same patient. One sample is from a kidney tumor biopsy. The other sample is a biopsy from the patient's other kidney, which seems to be perfectly healthy tissue.

Theoretically, you would expect that the two samples will have different amounts of certain messenger RNA transcripts, because the cells within the two tissue samples are performing different functions and are thus transcribing different genes. **You're curious to see which transcripts from the tumor are being produced at a significantly higher or lower number in the tumor tissue than the healthy tissue. These differences in transcript production in aggregate comprise the differential expression that you're studying.**

For modern differential expression studies, you're going to make sequencing libraries from the messenger RNA transcripts produced by the cells in each sample. You'll sequence those transcripts and then (after some cleanup of the raw sequencer data) align your sequenced reads to the genome. You'll see which regions of the genome have coverage, and compare the depth of coverage for the regions between samples to determine which regions have significantly different coverage depth. Then you can compare the transcripts of interest against a database of known transcripts to determine the genes (and the isoforms of the genes) that are being differentially expressed.

Now, in a real study you're going to ideally have multiple patients and tissue samples from several sites on each patient. You're not necessarily going to be studying tumors, either. There are lots of potential transcript expression patterns to be gleaned for pretty much any biological phenomenon through careful study design.
Question: Hi guys,

we are working on an university project where we want to find discriminating genes of different cancer types. For that we are using gene expression data of the TCGA dataset.

A naive approach would be to simply run some feature selection on tumor data of each type. However, we assume that this would not identify genes relevant for the tumor but the cell type itself. For example, we want to compare thyroid and lung cancer. Using only tumor data, we would expect that we find differentially expressed genes that are specific not for the tumor but for the original cell type itself. So we want to "normalize" thyroid tumor data with healthy thyroid tissue to find discriminating genes for thyroid tumor first that can now be compared with "normalized" genes for lung cancer.

We have some ideas how to do this ourselves but we suppose that this is not an uncommon task, so **has anyone heard of this "normalization" approach and how it usually is done**? We suppose that this needs to be done when clustering cancer types to see meaningful differences but we could not find this in the literature we read.

We hope we could state our problem in a comprehensible way, if not, feel free to ask. Thanks for your help!

Answer: Hey, you could try:

1. Obtain healthy / normal tissue specific expression data from one or more online databases ( see
here: https://www.biostars.org/p/90411/ - also look up *FANTOM5*)
2. Using these databases, determine which genes are specific to your
tissues of interest (e.g. thyroid). To do this, just do something
like converting the downloaded data (preferably it has a normal distribution) to Z-scores. Thus,
'tissue-specific' genes will have high Z-scores in the
tissues in which they are most expressed; whereas, ubiquitously
expressed genes will have low Z-scores

After that, when you conduct your differential expression comparison in tumours between thyroid and lung, you can just filter the list of differentially expressed genes for the genes encountered via the above 2-step approach) - this is the simplistic approach.

There should be a way, however, to actually 'normalise' your tumour expression data based on the results that you obtain from the tissue specific databases. For example, one could supply the healthy / normal tissue-specific data as priors to an empirical Bayesian regression model and then adjust your expression data based on these priors.

Also keep in mind that 'thyroid' and 'lung' refer to many different cell- and tissue-types.

Kevin
Question: Hi all,

I've written one perl script, which compares two vcf files and write output in third file. it is running good with small files, but when i am inputting big files in Gbs, my system becomes slow and it would hang.

I'm using dbSNP vcf file as input which is around 9GB.

Please Suggest me something that will make my perl script run faster. This is my first perl script and i'm new to perl.

Any help would be appreciated !!

Thank you !!!!

```
#!/usr/bin/env perl
use strict;
use warnings;
use Data::Dumper;
use List::MoreUtils;

open (FILE, "<", "dbSNP_in.vcf") or die "failed to open: $!\n";
my @array=();
my @CHR;
my @location;
my @rs;
my @ref_n;
my @alt_n;
foreach (@array)
{
chomp;
my ($chrom, $pos, $id, $ref, $alt, $qual, $filter, $info) = split(/\t/, $_);

push @CHR, $chrom;
push @location, $pos;
push @rs, $id;
push @ref_n, $ref;
push @alt_n, $alt;
}

open (FILE1, "<", "trial_rs.vcf") or die "failed to open: $!\n";
my @array1=();

open (OUT,">trial_output.vcf");
my @columns;

foreach (@array1)
{
chomp;
@columns=split(/\t/, $_);
my $i;
for ($i=0; $i<@array; $i++)
{
if (($columns[0] eq $CHR[$i]) and ($columns[1] eq $location[$i]) and ($columns[3] eq $ref_n[$i]) and ($columns[4] eq $alt_n[$i]))
{
$columns[2]=$rs[$i];
}

}

print OUT join("\t", @columns), "\n";
}
```

Answer: The problem is this:

my @array=(); ## slurp a whole Terabyte into RAM

Do not read a file like this, unless it is very small or your memory is huge, instead use:

```
while () { # read each line, then forget about it
chomp;
split /t/;
....
}
```

In addition your code is extremely inefficient:

You are iterating over all entries of a huge array for each line of a huge file. Hashes in Perl provide constant time access to elements given the hash key. Use a hash of hashes to store your filter table.

You seem to filter a large file by a small file, therefore:

- process the small file first (using while construct), parse the small file into a hash using the location as key of a nested hash
- process the large file as above and look up each of the entries in the hash created before

As a result you will only need as much memory as is required to store the small file.
Question: Dear all,

I have aligned a file with the command:

bwa mem -t 10 -o .sam

then I converted the SAM to BAM with:

samtools view -Sb .sam > .bam

and I wanted to make the index of it with:

samtools index .bam

but I got the error:

> [E::hts_idx_push] NO_COOR reads not in a single block at the end 0 -1
> samtools index: failed to create index for "output.bam"

To check if the error was due to missing alignment I checked the number of mapped reads with:

samtools -c -F 4 .bam

but I got the error:

> [main] unrecognized command '-c'

and when removing the flag -c I got:

samtools -F 4 .bam

> [main] unrecognized command '-F'

Why I got the [E::hts_idx_push] error? maybe a missing heading?

And why the flags are not recognized?

Thank you.

Answer: There are two things:

1) Indexing requires coordinate-sorted files.

bwa mem idx file_1.fq.gz file_2.fq.gz | samtools sort -o file_sorted.bam -
samtools index file_sorted.bam

Directly pipe the sam file to the sort, there is no need to keep the uncompressed and unsorted sam files. They only take up space.

2) Your -c -f -F etc... commands must be used like samtools view -f -F -c (...). You did not specify the subcommand they belong to.
Question: I have imputed genotype data from the Haplotyle Reference Consortium in VCF format, and would like to convert to BGEN 1.3 (or 1.2). I know that `qctool` version 2 support this, as does `plink` version 2.0.

I would like to inquire as to experiences anyone has had converting genotype probabilities (GP field in VCF) using these tools:

1. Is one tool preferred over the other?
2. Any caveats we should be aware of?

I would tend to use `plink`, as it is quite popular and actively developed, although the developers of `qctool` also developed the BGEN format, which may make it more ideal. I would assume these tools are equivalent for the task of converting VCF to BGEN, but your experience with this may suggest otherwise.

Answer: plink2 supports allelic dosages, not genotype probabilities. So if your VCF GP field says P(0/0)=0.2, P(0/1)=0.52, P(1/1)=0.28, plink2 will just save "dosage(1)=1.08", and then when it's asked to export to BGEN it'll encode the dosage as P(0/0)=0, P(0/1)=0.92, P(1/1)=0.08.

When this is fine, plink2 is a good choice. When this isn't good enough, use something like qctool.
Question: I have a list of gene symbols.

APOB, TTC39B, ATF3, RGS1, LIPG,

I am trying to get the genomic coordinates (as in bp) with +/-5KB position of these genes via UCSC TableBrowser / MySQL server using NCBI 36/hg18 build. I have tried to get this information via TableBrowser, but I can't find the chromStart and chromEnd field in the given table. Am I missing something ?

Also, please share your favorite tutorial / docs that explain the schema/tables in UCSC MySQL server.

Thanks in advance.

Answer: From [the table browser][1], select group= Genes , track= UCSC gene , table=knownGene and then '**describe table schema**'

You'll see that `knownGene` is linked to `kgXref`:

hg18.kgXref.kgID (via knownGene.name)

and `kgXref` contains a column named `geneSymbol`.

All in one, you can get the positions of the transcripts for those genes:

mysql -h genome-mysql.cse.ucsc.edu -A -u genome -D hg18 -e 'select distinct X.geneSymbol,G.chrom,G.txStart-5000,G.txEnd+5000 from knownGene as G, kgXref as X where X.geneSymbol in ("APOB", "TTC39B", "ATF3", "RGS1", "LIPG") and X.kgId=G.name'
+------------+-------+----------------+--------------+
| geneSymbol | chrom | G.txStart-5000 | G.txEnd+5000 |
+------------+-------+----------------+--------------+
| APOB | chr2 | 21072805 | 21125450 |
| ATF3 | chr1 | 210843616 | 210865739 |
| ATF3 | chr1 | 210800319 | 210865739 |
| ATF3 | chr1 | 210843616 | 210865704 |
| ATF3 | chr1 | 210849982 | 210864212 |
| LIPG | chr18 | 45337424 | 45378276 |
| LIPG | chr18 | 45337424 | 45367217 |
| RGS1 | chr1 | 190806479 | 190820782 |
| TTC39B | chr9 | 15156560 | 15302244 |
| TTC39B | chr9 | 15156560 | 15227442 |
| TTC39B | chr9 | 15171584 | 15302244 |
| TTC39B | chr9 | 15172968 | 15268702 |
| TTC39B | chr9 | 15172968 | 15227442 |
+------------+-------+----------------+--------------+

I agree that the information is hard to find: I only knew where to find it because I use to play regularly with those tables.

The **UCSC mailing list** is a good place to find this kind of information.

I also did a lot of reverse engineering by just 'greping' the flat files available from the UCSC.

[1]: http://genome.ucsc.edu/cgi-bin/hgTables?command=start
Question: Hi Biostars,

Can anyone helps me on the following.
I want to plot the expression values of miRNAs as function of samples using ggplot library.
I've a dataframe "reduced_data" containing 50 samples and 3 miRNAs + one column "class" (class of each patient: 6 classes).

reduced_data (ordered by class):

miRNA1 miRNA2 miRNA3 class
sample1 0.126766967 0.1396135 0.126766967 1
sample2 0.002176183 0.2699439 0.005832450 1
sample3 0.004802792 0.3232242 0.002334003 1
sample4 0.005626610 0.2865383 0.003429664 1
sample5 0.078054375 0.1274370 0.078054375 1
sample6 0.002946872 0.2648828 0.002946872 1
sample7 0.048704970 0.1786870 0.048704970 1
sample8 0.005917859 0.3780041 0.005917859 2
sample9 0.003999677 0.2946423 0.005693926 2
...

I want to use ggplot to plot the value of each miRNA as function of patients and classes, each miRNAx in one colour.

What I want to obtain (samples are stratified by classes):

![enter image description here][1]

What I've already tried:

mir1 <- data.frame(samples = rownames(reduced_data), class = reduced_data[,4], value = reduced_data[,1])
mir2 <- data.frame(samples = rownames(reduced_data), class = reduced_data[,4], value = reduced_data[,2])
mir3 <- data.frame(samples = rownames(reduced_data), class = reduced_data[,4], value = reduced_data[,3])
mirs <- list(miRNA1=mir1, miRNA2=mir2, miRNA3=mir3)
df_mirs <- cbind(cat=rep(names(mirs),sapply(mirs,nrow)),do.call(rbind,mirs))
ggplot(df_mirs, aes(samples, value, color=cat)) + geom_line()
# i got this error:
#geom_path: Each group consists of only one observation. Do you need to adjust the group aesthetic?

#By replacing samples by class in aes() function,
ggplot(df_mirs, aes(class, value, color=cat)) + geom_line()

I got the following:

![enter image description here][2]

Thanks in advance
[1]: http://www.student.montefiore.ulg.ac.be/~debit/ggplot/ggplot.png
[2]: http://www.student.montefiore.ulg.ac.be/~debit/ggplot/ggplot2.png

Answer: Then add a column sample to "reduced_data"
then

library("reshape2")
m <- melt(reduced_data, id.vars=c("samples","class"), variable.name=("miRNA"))
ggplot(m, aes(x=samples,y=value,colour=miRNA))+ facet_grid(. ~ class)

In general in ggplot2 its always bringing the data into the correct form:
a long table with one measured value per item
http://vita.had.co.nz/papers/tidy-data.pdf

And make your questions more reproducible
Question: I have RNA-seq data from two treatment groups (three replicates each) and I am trying to replicate the likelihood ratio test outlined on page 1516 of this paper by Marioni et al: http://www.ncbi.nlm.nih.gov/pubmed/18550803. So far I have been able to fit a Poisson regression model for each gene using the following function and command:

poissonGLM <- function(gene) {
data <- data.frame(treatment=pdata$treatment, count=gene)
model <- glm(count ~ treatment, data=data, family=poisson)
mean <- predict(model.pois, type="response"),
stdev <- predict(model.pois, type="response", se.fit=T)$se.fit
}

results <- apply(edata, 1, poissonGLM)

For each gene, this computes the expected mean and expected variance for each treatment group. However, I'm unsure how to use this model to compute the maximum likelihood estimates both under the null and alternative hypothesis (as described in the paper), and generate a P-value from the Chi-squared distribution.

Answer: anova(model, test="Chi") should give you the likelihood ratio test result comparing the effect of treatment vs intercept only model. This is the same as if you had manually fit a glm with treatment and a second glm with intercept only and computed their deviance and compared it to a chi-squared distribution (d.f. = 1 in your case).
Question: I'm working on a project where the read data consists of 6-mers ligated to the 5' end of a sequence. These 6-mers are important, but must be excluded during the alignment of the reads to which they are attached.

My current method of evaluating these experimental data consists of trimming off the 6-mers, aligning the trimmed reads, and then using the read query-name to go back and retrieve the 6-mer from the original FASTQ data.

This method currently works just fine. However, it would be far more efficient if I could use an aligner with the ability to soft-clip (not trim!) the reads, such that the clipped bases from each read aren't considered during the alignment, but will show up in the BAM data with the appropriate CIGAR soft-clipping designation.

Does anyone know if there's any aligner out there which supports this option? I've searched through several (bwa, bowtie, SOAP), but unless I'm misunderstanding the documentation, none of them support what I'm trying to do.

Answer: I would write a small script to move the 6bp to the FASTQ header:

```
my ($name, $seq) = ('', '');
while (<>) {
chomp;
if ($.%4 == 1) {
$name = @{split()}[0]; # save the read name
} elsif ($.%4 == 2) {
$seq = $_; # save the full sequence
} elsif ($.%4 == 0) {
# print the FASTQ header with the first 6bp in b6/q6 tags
print "$name b6:Z:" . substr($seq, 0, 6), "\tq6:Z:" . substr($_, 0, 6) . "\n";
# print sequence and quality with first 6bp trimmed
print substr($seq, 6), "\n+\n";
print substr($_, 6), "\n";
}
}
```

and then:

zcat myreads.fq.gz | perl move6.pl | bwa mem -C ref.fa - | gzip -1 > out.sam.gz

The output will have "b6" and "q6" tags. In general, option "-C" copies the FASTQ comment to the SAM output. You can use this approach to attach any meta information about reads in the SAM output. Bwa-aln also has a barcode option.
Question: Hey all,

I am following a protocol from a paper that uses the following pre-processing procedure:

> a. Read counts were normalized between samples using TMM (Robinson, M. D. & Oshlack, A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biology 11, R25 (2010)).

>b. Expression values for each gene were inverse normal transformed.

I used `edgeR::calcNormFactors` to normalize library size via TMM for part a, but I am confused on how to apply an inverse normal transform on my read counts together with my normalized library sizes. What is my misunderstanding? I know that I can apply other transforms like `cpm`, `rpkm`, etc., to the results of `calcNormFactors`, and it will transform using the normalized library sizes -- is there a similar function for inverse normal transformation?

Appreciate any help.

Answer: I am not convinced that the inverse normal transformation is a good idea, so we don't provide a function for it in edgeR. Here is how you would do it however:

```
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE)
n <- ncol(logCPM)
zvalues <- qnorm(ppoints(n))
z <- logCPM
for (i in 1:nrow(z)) z[i,] <- zvalues[order(order(z[i,]))]
```

The inverse normal values are now stored in `z`.