cannin/biostars_qa huggingface dataset:
Question: Hi all-

Apologies if this question has been asked before, but I couldn't find anything on biostars or elsewhere which quite captured what I am trying to do.

I have two sets of ChIP-seq peaks (WT and KO) that I have assigned to the nearest genes. So, I now have three lists of genes: one list for the genes which have a peak in both datasets, one list of genes that only appears in the WT dataset, and one list of genes that only appears in the KO dataset. What I would like to do now is perform some kind of differential gene ontology analysis between the peaks that appear only in the WT or only in the KO. In other words, my desired output would look something like "genes with peaks associated only in the WT sample are associated with GO term X, while genes associated with peaks only in the KO sample are associated with GO term Y, indicating that the KO sample loses peaks associated with term X and gains peaks associated with term Y." Can someone point me in the right direction?

Thanks!

Answer: You could give clusterProfiler a shot using `compareCluster()`. You'd need to prepare a data frame where each column is a gene list and then do something like:

require(clusterProfiler)
require(enrichplot)
require(org.Hs.eg.db)

cluster <- compareCluster(geneClusters = df, fun = "enrichGO", keyType = "YOUR_GENE_IDENTIFIER_TYPE", org = org.Hs.eg.db)

dotplot(cluster)
You'd need to swap out org.Hs.eg.db for another organism database if not human.

Vignette for clusterProfiler here: https://yulab-smu.github.io/clusterProfiler-book/chapter11.html
Question: Hello All,

This is a duplicate of a question I have posted on stack overflow but I think perhaps this community is better informed about possible tools that can solve it.

Question on stack overflow: http://stackoverflow.com/questions/32312790/converting-lists-of-data-into-informative-table-or-heatmap

So the issue is:

I am looking for a good way to visualise overlaps of genes found in multiple tests. I would like to check for overlap between multiple files containing lists of genes and output a table that shows which files contain which genes in a specific way (example output further below).

I have multiple text files with lists of genes. One gene per line. Files range from approximately 30-100 rows. To be as clear as possible I will show 4 example files that I have shortened for space.

File1:

```
NRG3
FOXP3
SHH2
ROBO1
PPP3CA
```

File2:

```
NRG3
SHH2
```

File3:

```
NRG3
ROBO1
```

File4:

ROBO1

I would like a way to take these files and create an output table that takes all the genes that are in the files and prints them as the rows of the first column (sorted alphabetically with each gene appearing only once). Then the following columns each represent an input file. If a gene appears in a file there will be an 'x' (or some arbitrary marker) to represent this in the relevant column. This would provide an easy way to visualise which genes appear in multiple files. Like this:

```
File1 File2 File3 File4
FOXP3 X
NRG3 X X X
PPP3CA X
ROBO1 X X X
SHH2 X X
```

It would be even more useful if, instead of an 'x' to represent if a gene appears in a file, this was shown in a heatmap color-gradient style way, so when a gene is found in only one file the relevant cell is shaded a light yellow, whereas if it appears in all files the cells are shaded a dark red. Are there R packages that exist to do this? However this is not essential, I just think it would be cool and improve the clarity of the visualisation.

I would appreciate any advice on how to go about doing this, especially if there are existing packages in R I am unaware about that do this already. Let me know how I can be more clear in explaining this problem.

Thank you for your help

Answer: This is how I would do it in R:

```r
FILES <- c("FILE1.txt", "FILE2.txt", "FILE3.txt", "FILE4.txt")
genes.all <- unlist(lapply(FILES, function(x){ readLines(x) }))
genes.all <- unique(genes.all)

res.df <- as.data.frame(lapply(FILES, function(x){
genes.all %in% readLines(x) }), row.names = genes.all)

colnames(res.df) <- FILES
```
Question: I have a list of about 1000 cancer cell lines with associated RNA-seq expression profiles per gene. I would like to be able to confirm that these cell lines have been verified or meet a standard that can be used in research since some cell lines can be contaminated or labeled incorrectly at some point. Is there a list of "verified" cancer cell lines somewhere? Or any advice on how to go about this?

Answer: There is a [Register of Misidentified Cell Lines][1]:

> ICLAC curates a Register of cell lines that are known to be
> misidentified through cross-contamination or other mechanisms (e.g.,
> mislabelling).

It seems to be actively maintained as well (the last update was less than a month ago).

[1]: http://iclac.org/databases/cross-contaminations/
Question: Hi all,

I've read that CNV-seq needs a reference hit file. Is this reference file the actual reference genome (the one I used to map my reads)? Or I have to pick up one of my samples to be a reference hit file ? (I don't have "treatment/control" samples). In the case that should be the actual reference genome, how do I get the .bam file? I mean, it's the reference genome, it doesn't make sense to map it against itself.

Thanks!

Answer: So far my solution was to create a "dummy" reference sample hit that contains all positions as they also do in CNVkit when no "normal" sample is provided.
Question: How can I best plot a histogram for billions of genotype quality values?

I have a simple one column file with billions of genotype quality values. The file is several GB uncompressed.

Is there a statistics library in Python or R that can build up a histogram by streaming trough the data?
Instead of loading everything to memory and then creating the histogram? I prefer using all of the data versus sampling it.

Or do I have to write a script first to collect the counts per bin and then give those count per bin to R for plotting? This functionality feels like it should already exist in a stats library some where.

I know the min and max of the values and would be able to specify a bin size.

Answer: *Or do I have to write a script first to collect the counts per bin and then give those count per bin to R for plotting? This functionality feels like it should already exist in a stats library some where.*

This is what I'd do. You're right that it might exist, but writing a 10 line perl script to bin them might be quicker than searching. I'd wager you'd be done before someone chimes in with the library to do it :)
Question: Hey everyone,

I am trying to transpose large datasets in excel but it keeps giving me the message we can't paste because copy area and paste area aren't the same size. Is there a way to transpose all the data at one time instead of piece by piece? One dataset has a great amount of rows and columns. Thanks.

Answer: Excel is not for Bioinformatics.

Use R to transpose the data matrix.

data <- read.table("your_file_location", sep = "\t", comment.char = "", stringAsFactors = F, header = T)

transpose_data <- t(data)
Question: Hello,everyone

I used Picardtools/MarkDuplicates as following:

java -jar picard.jar \
INPUT=A01-2m-D6.sort.bam \
OUTPUT=A01-2m-D6.bam \
METRICS_FILE=A01-2m-D6.metrics CREATE_INDEX=true 2>A01-2m-D6.picard_MarkDuplicates.stderr

then I met thus error:

picard.sam.markduplicates.MarkDuplicates done. Elapsed time: 2.38 minutes.
Runtime.totalMemory()=2058354688
To get help, see http://broadinstitute.github.io/picard/index.html#GettingHelp
Exception in thread "main" htsjdk.samtools.SAMFormatException: SAM validation error: ERROR: Record 78004, Read name E00477:238:hggvjccxy:8:2224:1519:63050, Zero-length read without FZ, CS or CQ tag
at htsjdk.samtools.SAMUtils.processValidationErrors(SAMUtils.java:454)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.advance(BAMFileReader.java:812)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:797)
at htsjdk.samtools.BAMFileReader$BAMFileIterator.next(BAMFileReader.java:765)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:576)
at htsjdk.samtools.SamReader$AssertingIterator.next(SamReader.java:548)
at picard.sam.markduplicates.MarkDuplicates.buildSortedReadEndLists(MarkDuplicates.java:495)
at picard.sam.markduplicates.MarkDuplicates.doWork(MarkDuplicates.java:232)
at picard.cmdline.CommandLineProgram.instanceMain(CommandLineProgram.java:269)
at picard.cmdline.PicardCommandLine.instanceMain(PicardCommandLine.java:98)
at picard.cmdline.PicardCommandLine.main(PicardCommandLine.java:108)

How can I deal with this problem ?Thanks

Answer: use https://broadinstitute.github.io/picard/command-line-overview.html#StandardOptions

VALIDATION_STRINGENCY=LENIENT
Question: I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

```
>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence
```

And I have a list like the following:

```
Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11
```

I want to retrieve those sequences that match with names present in the list.

Answer: **BASH**

```
IFS=$'\n'
for i in $(cat headersFILE)
do
line=$(grep -nr "$i" multiFastaFILE)
if [[ ! -z $line ]]
then
for j in $line
do
lineNumber=$(echo $j | cut -d':' -f1)
sed -n "$lineNumber p" multiFastaFILE
awk -v nb=$lineNumber 'NR > nb {if ($0 ~ ">") exit; else print $0 }' multiFastaFILE
done
fi
done
```

Where headersFILE is the file containing one header per line ; multiFastaFILE is obviously your multifasta file.

Every time a fasta sequence header of your multiFastaFILE contains a name described in your headersFILE, the header and its sequence will be displayed.
Question: Hi guys, I need to select a core of genes from an RNA-Seq dataset of 200 samples whose expression is able to explain an x-percent of variability of the entire dataset (i.e. across samples).
My point is that x-percent is normally an arbitrary threshold. Could be 75% or 80% according to what the scientist considers to be relevant. Is there a way to set in a quantitative way that threshold?
Through simulations for example? In other word is there a way to prove that threshold for example below 75% are no more informative to explain the variability of data?
In my mind I consider 75% a reasonable cut-off but I need to explain to non-statisticians and with numbers that values below are not of interest.

Thank you in advance

Answer: This can be done via PCA. Basically, if you perform PCA and then have 70% variation explained by PC1, you can then extract the variables / genes that are responsible for the majority of this [70%] variation along PC1. If your PC1 explains a lower number, like 25%, then you may have to look at numerous PCs that account for, e.g., up to 80% variation.

Curiously, I have implemented functionality for this in *PCAtools* (released with Bioc 3.9), and Aaron Lun has also recently added functionality for identifying the ideal number of PCs to retain.

For Aaron's code, look here: Determine optimum number of PCs to retain

For identifying variables responsible for variation along each PC, look here: Determine the variables that drive variation among each PC

This Bioc package is new, so, please make any recommendations for its future development.

Kevin
Question: Hi all

I called six SNP's files as individual, I want to merge them such that considering the position and location. I want to do that for converting them as integer numbers 0,1,2.

The question is:

Could anyone please help me how I can merge them as following?

REF is hg19 , ALT1 is first patient, ALT2 second patient ... so on till ALT6 sixth patient.

```
#CHROM POS REF ALT1 ALT2 ALT3 ALT4 ALT5 ALT6
chrM 3 T C G A C T C
chrM 4 C A C T A G C
chrM 150 T C T C C G A
chrM 195 C T C T C A T
chrM 410 A T T C C T C
chrM 711 G A C T T G G
chrM 1890 G . C T C A C
chrM 2354 C T T C A G C
chrM 2485 C T A G G A C,T
chrM 3457 T C G A G A C
chrM 4162 C T T A T C,A A
chrM 4217 T C G T A G T
chrM 4918 A G C . G A A
chrM 5581 C T G A A G .
chrM 8698 G A G A A C A
chrM 8702 G A G C G C A
chrM 9378 G A C T G A C
chrM 9541 C T C T C T C
chrM 10284 A G G A A C C
chrM 10399 G A G A A G T
chrM 10464 T C C G T C G
chrM 10820 G A G T . C A
chrM 10874 C T G T G C,T G
chrM 11018 C T C T A C C
chrM 11252 A G . C G A T
chrM 11723 C T . A C T T
chrM 11813 A G G A C A C
```

Is that possible? I wrote period because someone told me you should have these periods if the positions there!

Thank you in advance

Bing

Answer: Assuming you have 6 VCF files, you can use GATK.

java -jar GenomeAnalysisTK.jar \
-T CombineVariants \
-R reference.fasta \
--variant input1.vcf \
--variant input2.vcf \
--variant input3.vcf \
--variant input4.vcf \
--variant input5.vcf \
--variant input6.vcf \
-o output.vcf \
-genotypeMergeOptions UNIQUIFY

Then you can use `VariantsToTable` to turn output.vcf into a table as requested:

java -jar GenomeAnalysisTK.jar \
-R reference.fasta -T VariantsToTable \
-V output.vcf \
-F CHROM -F POS -F ID -F QUAL -F -GF GT \
-o output.table
Question: Hi,

I have been trying to download ChIP-seq data using SRA-toolkit.

prefetch -v SRR4032857

I get an error

> Maximum file size download limit is 20,971,520KB
>
> 2017-01-03T09:32:45 prefetch.2.3.5 err: error unexpected while
> resolving tree within virtual file system module - failed to resolve
> accession 'SRR925811' - Obsolete software. See
> https://github.com/ncbi/sra-tools/wiki ( 406 )
>
> Redirected!!!
>
> 2017-01-03T09:32:46 prefetch.2.3.5 err: name incorrect while
> evaluating path within network system module - Scheme is 'https'
> 2017-01-03T09:32:46 prefetch.2.3.5 err: path not found while resolving
> tree within virtual file system module - 'SRR925811' cannot be found.

I checked my hard disk space and also changed the default download folder from ./vdb-config -i. I downloaded the SRA-toolkit version 2.8.1.

What could be the issue and how can i resolve it?

Answer: As per the link in the error message, SRA has moved from HTTP to HTTPS, rendering all previous versions of the SRA toolkit obsolete. You say you downloaded 2.8.1, but you still appear to be running 2.3.5, so I'd check what version is in your `$PATH` and where you installed the latest version.
Question: Hello,

I am using PLINK2 to create a profile score and am confused to what the output represent:

My score file has ~84k SNPs and my profile output looks like this

```
FID IID PHENO CNT CNT2 SCORE
00 1 1 173748 83462 0.017818
01 1 1 173740 83216 0.017095
02 2 1 173728 83738 0.017375
```

Plink says

```
CNT -Number of non-missing SNPs used for scoring
CNT2 -The number of named alleles
```

It seems to me that plink might have interchanged CNT and CNT2 column but I do not understand what "number of named alleles" mean ? Also, the scoring is such, all my individuals are scored in the range of 0.017095 to 0.017999, isn't that a bit strange ?

Many thanks

Answer: Your score file should have lines looking like

rs4970405 G 0.0012

The second column names an allele which contributes to the score; CNT2 is the total number of named alleles actually observed. CNT is the number of *alleles* checked (the PLINK 1.07 documentation is incorrect here), so it's always at least as large as CNT2.

The small range is unsurprising if most of your scored variants have low minor allele frequencies. You can use the 'center' modifier to normalize them to have mean zero.
Question: I have RNAseq reads (single-end, 49-bp), I would like to do differential expression analysis. I know the standard pipeline is some rnaseq aligners ( e.g. Tophat, soap etc,) followed by featureCounts/Cufflinks(with assembly)/HTSeq and EdgeR/DESeq.

My question is,
1) will it be appropriate to use BWA aligner instead of other standard RNAseq specific aligners? And
2) why there is lot of difference in percentage of mapping reads when you do alignment with BWA and other RNAseq specific aligners?

Answer: BWA isn't going to handle spliced alignment terribly well, so it's not normally recommended for RNAseq datasets. That's also why it gives lower alignment rates.

BTW, give STAR a try. It's much faster than tophat.
Question: Hi everyone,

Sometimes I am asked to mention the most exciting moment or part of my research works, I often say that I was excited first time I know how to work with command lines in Linux. However I think might be there is a specific purpose behind asking this question. What would you reply if you asked so?

Thank you very much

Answer: I'd say fixing a bug/error in a script you wrote, or significantly speeding things up is particularly rewarding. But also getting nice feedback from users is great.
Question: Hi,

Does anybody know of any script/tool which could convert a MAF file to VCF?

https://wiki.nci.nih.gov/display/TCGA/Mutation+Annotation+Format+%28MAF%29+Specification?showChildren=false

Thanks,
Tinu

Answer: There's a script called maf2vcf in the vcf2maf repository.
Question: I want to transform a TCGA mRNA expression matrix (in linear data format) to log2-ratios and then run a feature (gene) selection, selecting the 1000 most variant genes (genes with higher standard deviation across samples). The workflow is the following:

1. Select "good" genes before log2ratio (genes each with median signal at least t in p% of samples);
2. On selected genes, run log2ratio, dividing each gene by its median signal and then log2-transforming the result matrix;
3. Select the 1000 most variant genes along all samples.

How do I select t and p?

Answer: There is not a general solution to select "t" and "p". Such choices are largely arbitrary. Furthermore, for an array platform, if one assumes that "t" has something to do with "expressed", the value for "t" will differ for each probe on the array.

Since you are ultimately going to filter based on variance, I'd suggest starting with your median-centered, log-transformed data and simply choose the top 1000 most variable genes.
Question: Hi there,

I'm trying to create a pheatmap with annotationcolors based on groups. I'm an sociologist, so this might not be the common way to use pheatmap, but I found it to be great at displaying social groups as well. But,, when I try to plot an assymmetrical matrix, where the rows and columns do not match, I get an error:

```r
Error in check.length("fill") :
'gpar' element 'fill' must not be length 0
```

The issue is also tackled [here][1], but as I understand it, the solution is to make sure rownames and colnames are identical. Which isn't possible for my purpose.

here is an reproducible example:

```r
DAT <- matrix(c(
771,1176,1018,1042,275,949,823,1250,348,622,1244,
465,984,766,790,573,784,592,1250,1222,421,1103,
171,288,194,80,547,528,249,959,72,81,224,
132,122,90,50,44,206,127,108,52,59,144,
0,9,5,9,0,10,17,0,0,0,0),
nrow = 5, ncol = 11, byrow = TRUE)

group <- c("grp1","grp2","grp3","grp4","grp5","grp6","grp7","grp8","grp9","grp9","grp9")

annotation_c <- as.data.frame(group)

colnames(DAT) <- paste(group," ", "no.",seq_len(11),sep="")
rownames(DAT) <- paste(rep(c("grp"),5),91:95,sep="_")

# annotation_r <- as.data.frame(rownames(DAT))

pheatmap::pheatmap(DAT, color = colorRampPalette(rev(brewer.pal(n = 7, name = "RdYlBu")))(length(seq_len(10))), cluster_cols=FALSE,cluster_rows=FALSE ,
annotation_col=annotation_c)
```

Thank you for your time,
Emil

[1]: https://www.biostars.org/p/175858/

Answer: As E.Rempel wrote:

*I added rownames to the annotation_c data frame*

> rownames(annotation_c) <- colnames(DAT)

*and your script worked.*
Question: Hello all,

Recently, I started to work with alignment data.
At the moment I have alignments for Ribo and RNA seq data that I can visualize with IGV tool.
I am trying to understand the sam files of these alignments to then to try to get some information from them.

Currently, I have two questions:
The first one is to understand the values of MAPping Quality.

I observed this case:

In my reference genome I have this sequence from chromosome 1.

>chr1:10538-10568
CCACCGAAATCTGTGCAGAGGACAACGCAGC

Into the SAM file of the alignment of RiboSeq (single end) reads, I got this:

SRR4293695.162674600	272	chr1	10538	0	30M1S	*	0	0 CCACCGAAATCTGTGCAGAGGAGAACGCAGA FFFFAFFFFFFF7chr1:10538-10569
CCACCGAAATCTGTGCAGAGGACAACGCAGC --> Sequence in genome
CCACCGAAATCTGTGCAGAGGAGAACGCAGA --> Alignment sequence

Why this happened? Or why one difference is counted as Match?

Thank you in advance for your answers!

Answer: A mapping quality of zero usually indicates a multiply mapped read (one that aligns equally well in multiple locations). The alignments often indicate the other locations, but in this case, they do not. The flag 272 has more information

0x110	272	REVERSE,SECONDARY

See how this alignment is a secondary alignment.

The sequence is aligned with `30M1S` with 30 matches and one soft clip. Th last base is chopped off and is not part of the 30M aligned region. This works this way since the end clip penalty is smaller than a mismatch penalty. There are specific instances where using this type of scoring is advantageous.

That being said, the `30M` would also include mismatches if you actuall had them as the `M` stands for match or mismatch.

It is quite misleading actually, the mismatching alignments would be identified from the `MD` and `NM` tags of the alignment. But since both of these tags are optional it is very much possible to have an alignment where you can't easily tell by eye if the alignment has mismatches. Umm ... yeah ... I agree that's just plain ridiculous.

The `samtools calmd` command will fill in missing MD tags.