# Question
Title: How to select the most representative pathways from a gene enrichment analysis?
I have perform an enrichment analysis to a cluster of genes. The output is a list of pathways and their p-value (the pathways are selected because p-value \< 0.05). The list is still quite long, so I want to reduce it. For that purpose I have a calculated the Dice coefficient of the pathways in a matrix $p$x$p$ where $p$ is the number of pathways in the list. I want both the ones that are more different (they overlap less, their Dice coefficient is lower) and the pathways more representative of the most similar pathways (So if a there is a group of 5 pathways that overlap over 0.8 take just one).

How can I select the most representatives pathways?

There is a similar tool for GO but it relays on discarding not significant GO, while here all the initial pathways are already significant.

If I do a clustering of the pathways using the Dice coefficient matrix I don't know where (or how) to cut.

I tried using the height to select the pathways. But I am unsure of the interpretation of height.

Some other tools I have seen use a multidimensional scaling plot, but I am not sure if performing it and cutting at certain point of the first dimension would help.

# Answer
One way to group similar significant pathways, is to quantify how many genes overlap between the pathways, and then use this in clustering (heatmap). I have made a tool in R which calculates the overlap index between GO terms and subsequently clusters them in a heatmap. Overlap index is the fraction of genes that overlap (number between 0-1). Also Pearson correlation can be used for clustering instead of the overlap index. My package (gogadget) works only with goseq analysis, but you can use goseq also for reactome or kegg data.

I have used this overlap clustering approach for different data sets now, and usually we can reduce 200-300 GO terms into 10-20 functional groups.
> 3 votes

# Answer
This sounds like something that might be amenable to a clustered heatmap plot, or a correlation matrix plot, or something similar. Have you looked at a correlation matrix of the dice coefficient matrix (or maybe just a heatmap plot of that matrix without the correlation matrix)?

The `corrplot` package looks like it might be useful, in particular the `hclust` / drawing rectangles presentation.

I can't vouch for this package though; it's just something I found by a search for "R plot correlation matrix".
> 4 votes

# Answer
If you're happy with a more confident ranking of the most representative gene sets, rather than necessarily cutting down the list, you might try EGSEA. It uses an ensemble approach to give a ranking of the most relevant gene sets, and also produces an interactive HTML output with statistics, heatmaps, pathway maps, summary plots and GO graphs which allows you to examine the output at varying levels of granularity.

You can read the paper on bioRxiv or download the package from Bioconductor.
> 2 votes

---
Tags: pathway
---


# Question
Title: Is it possible to check if a patient has the HLA-B27 antigen from his exome stored in a VCF file
Is it possible to check if patient has the HLAB27 antigen by comparing his exome stored in VCF with all known HLAB27 variants?

For instance, rs4349859 is a known HLAB27 SNP.

If I can find rs4349859 in the patient's VCF file, could I say he is HLAB27+ ?

# Answer
> 2 votes
Some papers claim that you can do HLA typing from exome sequencing data. I have not validated this myself. I also guess it depends on capture. Nonetheless, I don't think you can type HLA with VCF generated with a standard pipeline, which is not specifically designed for HLA. Current variant callers won't work on HLA genes due to their high diversity. Also in general, there are thousands of HLA-B alleles. In most cases, you can't find the allele by looking at only one SNP.

No, you can't check if the patient has HLA-B27 only from VCF. There might be a hope if you get raw reads and run a proper HLA typing tool designed for such data.

---
Tags: vcf, exome
---


# Question
Title: Extracting features or gene from PCA after calculating PCA for downstream analysis
The quote below is from this paper:

> We performed principal component analysis (PCA) of low-coverage sequencing data to identify genes explaining variation across cells. PCA separated the cells into groups corresponding to the source populations (Fig. 2c and Supplementary Figs. 3–5), and the genes distinguishing each group reflected the biological properties of the cell types (Supplementary Fig. 5 and Supplementary Table 3). PCA of low- and high-coverage sequencing data revealed a remarkably similar graphical distribution of analyzed cells, and the majority (78%) of the top 500 genes determined by PCA were shared between PCA performed on low- and high-coverage data (Supplementary Figs. 4 and 6 and Supplementary Table 4).

How are they performing PCA and then finding out the components which contribute to the strongest PC and then taking out the genes for downstream analysis? Let's say there are genes involved in PC1 and PC2 which define a certain cell type and that distinguish it. I think for a biologist, like me, I just want to get the genes or gene list that are getting involved in determining lineage or cell type. How can I do this?

As a to test what they might be doing i did this

```
test<- read.csv('NON_CODING.csv',header = T,row.names = 1)

# preform PCA
pca = prcomp(t(test), center=TRUE, scale=TRUE)

```

then after calculating PCA i do see this in my console when I type pca

```
> pca
Standard deviations (1, .., p=16):
 [1] 3.826851e+01 2.080405e+01 1.568739e+01 1.349256e+01 1.119348e+01 9.980547e+00 8.365034e+00
 [8] 8.098841e+00 7.519507e+00 6.184505e+00 5.880260e+00 5.139609e+00 4.851091e+00 4.335606e+00
[15] 3.870918e+00 3.173450e-14

Rotation (n x k) = (2899 x 16):
                                       PC1           PC2           PC3           PC4
5S_rRNA                      -1.090574e-02 -2.412665e-03  1.637689e-02 -3.603865e-02
AB019441.29                  -1.928250e-02  1.821083e-03 -9.713724e-03 -1.978737e-03
ABBA01017803.1               -1.823266e-02 -3.727144e-03 -9.790131e-04 -3.937062e-02
ABC14-1080714F14.1            2.438024e-02  3.816019e-03  1.657784e-04 -7.141480e-03
ABC7-481722F1.1               2.467403e-02  2.873432e-03  3.781153e-03  3.790824e-03
AC000036.4                   -1.066777e-02  2.838305e-02  2.642673e-02 -7.998608e-03
AC000089.3                   -1.249026e-02 -7.870864e-04 -2.569602e-02 -3.165073e-03
AC000120.7                   -1.313696e-02  5.923245e-03  1.633255e-02  4.077602e-02
AC000123.4                    2.415996e-02  1.029732e-02  8.930792e-03 -1.118660e-02
AC000403.4                    8.692835e-03 -3.060170e-02  3.746839e-02 -6.968971e-03
AC002064.4                   -1.799590e-02  2.147160e-02 -2.303742e-02  7.647937e-03

```

Now how do i find which are the list or set of genes that gives me PC1 or PC2 so on

Any suggestion or help how to get the genes that make most difference between two major component that would be very helpful

# Answer
The help page explains that there is a component which:

> rotation
> the matrix of variable loadings (i.e., a matrix whose columns contain the eigenvectors). The function princomp returns this in the element loadings.

So using using the rotation one can use it to calculate the importance of each feature in the input:

```
C <- chol(S <- toeplitz(.9 ^ (0:31))) # Cov.matrix and its root
set.seed(17)
X <- matrix(rnorm(32000), 1000, 32)
Z <- X %*% C  ## ==>  cov(Z) ~=  C'C = S
pZ <- prcomp(Z)
genesImportance <- Z %*% pZ$rotation
genesImportance2 <- pZ$x
all.equal(genesImportance, genesImportance2)
## [1] "Mean relative difference: 0.03680457"

```

This new matrix has the loading for each gene for each component. So you can now select them according to how much do they participate in each component or whatever you want. Or you can use it to predict the position of samples you didn't set in the initial data.
> 7 votes

---
Tags: r, rna-seq
---


title: Get gene annotation type from gene name

text:
# Question
Title: Get gene annotation type from gene name
I am having `RNA-seq` data set for `C.Elegans`. I do not know which annotation scientists used for naming the genes of it, but I need to search somehow for the `^ERCC` and `^mt-` genes like specified in the tutorial:

https://www.bioconductor.org/help/workflows/simpleSingleCell/

If I just implement `grepl` search for the `character` array of gene names that I have, I am always getting all values `FALSE` as if there are no `spike-in`s and mitochindrial genes.

Any suggestion would be greatly appreciated.

Sample gene name from the dataset: `WBGene00023193`

# Answer
> 2 votes
The tutorial you linked to does not mention ERCC or mt genes.

Before running `grepl` or any other search, you should first check your inputs. It looks like you have WormBase gene IDs if all your genes look like `WBGene00023193`. Therefore, your question is actually how to convert WormBase gene IDs to gene symbols.

Here is one potential solution: https://github.com/WormBase/website/issues/4297

After you convert, you have to make sure that you have ERCC spike-ins in your gene table (probably not) and you mitochondrial genes start with "mt-".

---
Tags: rna-seq, gene
---


# Question
Title: Which tools are available to call somatic mutations on non-coding regions from whole genome sequencing data?
I have a set of whole genome sequencing samples, some of which have matched normal while some not. I want to call somatic mutations on non-coding regions.

I was looking for GATK best practices (like this) but it seems it is not yet ready for production. Which alternative tools do you suggest?

# Answer
> 1 votes
Those will have to be treated as separate analysis. You can use Mutect or VarScan for the matched samples and the tools listed in the other question for the unmatched samples however that means that you won't be able to compare the results easily. If the inter-sample comparison is important you will have to run all samples as unmatched. Personally I have used LoFreq for somatic variant detection with some success. It is very sensitive so the false positive rate can be high depending on how noisy your data is.

---
Tags: wgs, mutations, somatic, non-coding
---


# Question
Title: Mini-bulk RNA-seq and scRNA-seq analysis differences?
We have RNA-seq data from libraries prepared using the Smartseq2 single-cell protocol on **500 cells** (mini-bulk) / library. The complexity is much better than with single-cells (~14k genes for 1.5M reads).

6 cell types / biological replicates were collected with flow cytometry, and there are 4 control and 5 diseased biological replicates (patients).

We are mostly interested in the differences between control and disease either overall or within each cell type.

Can single-cell tools (e.g. Seurat) be used on this dataset? It is neither bulk nor single-cell data. How would you approach this?

# Answer
> 2 votes
Single-cell tools were developed to address some specific factors unique to scRNA-seq, such as:

* a lot of dropouts, so hundreds or thousands of detected genes, not tens of thousands (for mammalian genomes)
* hundreds or thousands of replicates
* not necessarily well-defined groups

Thus, I would say your experiment is a lot closer to bulk RNA-seq than scRNA-seq.

---
Tags: rna-seq, sequence-analysis, single-cell
---


# Question
Title: Where can I get the population allele frequency vcf file?
I want to use GATK to estimate cross-sample contamination for Whole Genome Sequencing data.

The specific tool is ContEst and it is run with:

```
java
 -jar GenomeAnalysisTK.jar \
 -T ContEst \
 -R reference.fasta \
 -I:eval tumor.bam \
 -I:genotype normal.bam \
 --popfile populationAlleleFrequencies.vcf \
 -L populationSites.interval_list
 [-L targets.interval_list] \
 -isr INTERSECTION \
 -o output.txt

```

The `--popfile` option requires a vcf file representing population allele frequencies. Where can I get it? Is it the one I can get from dbSNP?

I'm using the human\_g1k\_v37.fasta reference, a version of hg19, "for the phase1 analysis we mapped to GRCh37. Our fasta file which can be found here called human\_g1k\_v37.fasta.gz, it contains the autosomes, X, Y and MT but no haplotype sequence or EBV".

# Answer
> 5 votes
On the GATK forum they've recommended the population stratified VCF file for this purpose.

---
Tags: vcf, variant-calling, gatk, quality-control, wgs
---


# Question
Title: Downloading URL using AWK for fixed fields
I have a text file as follows

```
a   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz
b   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz
c   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz

```

The file is tab delimited and has 8 columns with 8th column being the url link I need to download the files one by one using wget and changing the name of the downloaded file with -O option in wget. The argument I want to give to the -O option is combination col1 and col2 separated by _ so that for first line I get the downloaded file as a\_k.fastq.gz and so on. For example, the first line of the file above should give:

```
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

```

Following the post from here I realized looping through the file is slow and sometimes difficult. How can I download these files using wget with awk or cut as the tool for fixed fields

# Answer
```
awk '{print "wget " $9 " -O " $1 "_" $2 ".fastq.gz"}' tmp

```

will give you

```
wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz -O b_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz -O c_k.fastq.gz

```

If you add a `| bash` to the command it will then execute all the `wget` commands.
> 5 votes

# Answer
A variation of Bioathlete’s answer that doesn't suffer from escaping issues:

```
awk 'BEGIN{ FS="\t" } { printf("%s\t%s_%s.fastq.gz\n", $8, $1, $2) }' input |
while IFS=$'\t' read -r url filename; do
    wget -O "$filename" -- "$url"
done

```

Or entirely with shell parsing:

```
while IFS=$'\t' read -r c1 c2 c3 c4 c5 c6 c7 url; do
    wget -O "${c1}_${c2}.fastq.gz" -- "$url"
done < input

```

Or entirely without shell parsing but with GNU `xargs`:

```
awk 'BEGIN{ FS="\t" } { printf("%s_%s.fastq.gz\n%s\n", $1, $2, $8) }' input |
xargs -r -d '\n' -n 2 -- wget -O

```

… you could even run multiple `wget` instances in parallel from GNU `xargs` with the `-P` option, though you should set the `-q` option on `wget` to disable progress bars overwriting each other on the same terminal. Then again it would be better to use a real command-line download manager like Aria or Axel for that task.
> 4 votes

---
Tags: awk, data-download
---


# Question
Title: Determine if a gene is mitochondrial or not
I need to determine whether a gene is mitochondrial or not for C. Elegans automatically from its name by running a regular expression on a dozen of thousands of gene names. Currently I am having gene symbols like:

```
clec-190, clec-189, 21ur-8912, nduo-1, 21ur-6539, 21ur-13096, nduo-6, 21ur-14241, ndfl-4

```

I tried running `grepl("^mt-", symbols)` and it found nothing. I am sure that there are some mitochondrial genes just from their descriptions. I would assume that first I need to convert gene symbols to some other symbols to do that but I am not sure here.

Basically, I am trying to do it in a similar way to the tutorial of single cell RNAseq.

# Answer
You can't determine if a gene is mitochondrial purely from its name. The simplest route would be to use the wormbase biomart and either download the gene names and their associated chromosomes or use the biomaRt package from bioconductor to do the same thing. You can upload a list of the genes you're interested in querying in either case (or just download all of them, and filter).

From your list, only `nduo-1`, `nduo-6`, and `ndfl-4` are mitochondrial.
> 10 votes

---
Tags: r, rna-seq
---


title: Read counts from BAM file

text:
# Question
Title: Read counts from BAM file
I have few BAM files which are generated from Ion Torrent Server (ampliseq) aligned to hg19 genome. I want to extract read counts from the bam files and I know that "featureCounts" can be used for this. But before getting the read counts, to view the bam file I gave the following command:

```
samtools view -H IonXpress_input.bam

```

Output:

```
@HD VN:1.4  SO:coordinate   GO:none
@SQ SN:NR_039978    LN:984
@SQ SN:NM_130786    LN:1766
@SQ SN:NM_138932    LN:9293
@SQ SN:NM_033110    LN:2678
@SQ SN:NM_000014    LN:4678
@SQ SN:NM_144670    LN:5192
@SQ SN:NR_040112    LN:1201
@SQ SN:NM_017436    LN:2106
@SQ SN:NM_016161    LN:1771
@SQ SN:NM_015665    LN:1854

```

I guess these are from header. After these I see something like following:

```
@RG ID:JODXM.IonXpress_022.1    CN:S5TorrentServerVM/MT205524   DT:2017-12-07T12:04:52+0100 FO:TACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGATCGATGTACAGCTACGTACGTCTGAGCATCGA KS:TCAGTTCGAGACGCGAT    PG:tmap PL:IONTORRENT   PU:s5/540/IonXpress_022 SM:SNU449_TEAD4

```

Haven't seen a bam file like this before. Any idea what are these and how I can get read counts?

# Answer
Your data is not aligned to hg19, but to a bunch of RNA ref sequences. If you would align to hg19 you'll get each chromosome instead of the NR\_\* or NM\_\* accession codes with your `samtools view -H` code. With `samtools idxstats file.bam` you'll get the reads per chromosome (or per nucleotide sequence in your case).

For use in `featureCounts` my advice is to get the fastq files, and align it yourself to hg19 (or hg38 to be more up to date).

If your data is properly aligned to hg19 your output should be something like this:

```
samtools view -H aligned2hg19.bam

@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:249250621
@SQ SN:chr2 LN:243199373
@SQ SN:chr3 LN:198022430
@SQ SN:chr4 LN:191154276
@SQ SN:chr5 LN:180915260
@SQ SN:chr6 LN:171115067
@SQ SN:chr7 LN:159138663
@SQ SN:chr8 LN:146364022
@SQ SN:chr9 LN:141213431
@SQ SN:chr10    LN:135534747
@SQ SN:chr11    LN:135006516
@SQ SN:chr12    LN:133851895
@SQ SN:chr13    LN:115169878
@SQ SN:chr14    LN:107349540
@SQ SN:chr15    LN:102531392
@SQ SN:chr16    LN:90354753
@SQ SN:chr17    LN:81195210
@SQ SN:chr18    LN:78077248
@SQ SN:chr19    LN:59128983
@SQ SN:chr20    LN:63025520
@SQ SN:chr21    LN:48129895
@SQ SN:chr22    LN:51304566
@SQ SN:chrX LN:155270560
@SQ SN:chrY LN:59373566

```
> 7 votes

---
Tags: rna-seq, sam, samtools
---


# Question
Title: Duplicate gene symbol handling in GEO gene expression data
I have downloaded a gene expression data from GEO database (GSE3268) in which in some of its rows there are duplicate gene symbols.

For example **TP53** exists in two rows with different expression values and different *GenBank Accession* IDs.

They have the following *Nucleotide Title*:

> Human p53 cellular tumor antigen mRNA, complete cds
>
> Homo sapiens tumor protein p53 (TP53), transcript variant 1, mRNA

I want to map these data with edge list of gene symbols to make a network. Which expression value should I consider for such duplicate gene symbols for mapping?

# Answer
Since this is a microarray (U133A), what you're seeing is that each probe is associated with one or more transcripts. The general strategy when dealing with microarrays is to use RMA (or fRMA) to summarize to probesets and then `collapseRows` to summarize to gene-level information. Since this is all typically done in R, you can find some further discussion on the bioconductor site. I would suggest going through this process and then only using the result of that in building your network (whether you're happy with sticking to transcript-level quantifications or gene-level quantifications for building your network is up to you).
> 5 votes

---
Tags: database, gene
---


title: sorting BAM file error using samtools

text:
# Question
Title: sorting BAM file error using samtools
I have few bam files and would like to get read counts using

```
samtools idxstats

```

\[Data is aligned to hg19 transcriptome\]. To use that command I need a sorted bam file. So to sort them I gave the following command.

```
samtools sort -T /tmp/input.sorted -o input.sorted.bam input.bam

```

This ended up showing:

```
[W::bam_hdr_read] EOF marker is absent. The input is probably truncated.
[E::bgzf_read] bgzf_read_block error -1 after 0 of 4 bytes
[bam_sort_core] truncated file. Aborting.

```

And the tmp dir has 6 input.sorted.nnnn.bam files.

# Answer
The cause of the error is that your input file is truncated or /tmp is running out of space. If you can do `samtools view -H input.bam` without error (reading the header also checks for the magic number at the end of the file) then it's like the that /tmp is filling up.
> 4 votes

---
Tags: rna-seq, sam, samtools
---


# Question
Title: Exception : [Failed to read 8 bytes from input stream! Read 0] salmon quant was invoked improperly
I am trying to run salmon and it keeps giving me 2 java exceptions:

```
java.io.EOFException
    at java.util.zip.GZIPInputStream.readUByte(GZIPInputStream.java:268)
    at java.util.zip.GZIPInputStream.readUShort(GZIPInputStream.java:258)
    at java.util.zip.GZIPInputStream.readHeader(GZIPInputStream.java:164)
    at java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:79)
    at java.util.zip.GZIPInputStream.<init>(GZIPInputStream.java:91)
    at uk.ac.babraham.FastQC.Utilities.MultiMemberGZIPInputStream.<init>(MultiMemberGZIPInputStream.java:37)

```

and

```
java.lang.ArrayIndexOutOfBoundsException: -1
    at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.calculateDistribution(SequenceLengthDistribution.java:100)
    at uk.ac.babraham.FastQC.Modules.SequenceLengthDistribution.raisesError(SequenceLengthDistribution.java:184)
    at uk.ac.babraham.FastQC.Report.HTMLReportArchive.startDocument(HTMLReportArchive.java:336)
    at uk.ac.babraham.FastQC.Report.HTMLReportArchive.<init>(HTMLReportArchive.java:84)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.analysisComplete(OfflineRunner.java:155)
    at uk.ac.babraham.FastQC.Analysis.AnalysisRunner.run(AnalysisRunner.java:110)
    at java.lang.Thread.run(Thread.java:745)

    at uk.ac.babraham.FastQC.Sequence.FastQFile.<init>(FastQFile.java:80)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:106)
    at uk.ac.babraham.FastQC.Sequence.SequenceFactory.getSequenceFile(SequenceFactory.java:62)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.processFile(OfflineRunner.java:129)
    at uk.ac.babraham.FastQC.Analysis.OfflineRunner.<init>(OfflineRunner.java:102)
    at uk.ac.babraham.FastQC.FastQCApplication.main(FastQCApplication.java:316)

```

I installed `salmon` using `Anaconda`s' `conda install -c bioconda salmon` and all other necessary packages in the same way. I thought that it could be a problem with `fastqc`, so I uninstalled it and then installed it manually (through the fastqc`.zip` file), but the output remained the same. I am using `Java8`:

```
java -version
openjdk version "1.8.0_92"
OpenJDK Runtime Environment (Zulu 8.15.0.1-linux64) (build 1.8.0_92-b15)
OpenJDK 64-Bit Server VM (Zulu 8.15.0.1-linux64) (build 25.92-b15, mixed mode)

```

I am running the script on a cluster with `SLURM`. Locally on my laptop the script is running without any issues. I am really stuck and do not know what to try, so any suggestions would be greatly appreciated.

`Update` Another issue that I just found is the following (the last command that `salmon` is trying to execute):

```
Exception : [Failed to read 8 bytes from input stream! Read 0]
salmon quant was invoked improperly.

```

Then I searched online for the issue and found the following thread:

https://github.com/COMBINE-lab/salmon/issues/129

However, updating `salmon` to the newest `0.9.1` version did not solve the issue.

`Edit` The command I am running is the following one:

```
python "$path/salmon_rna_seq.py" -L 'IU' -E "$path/jyvo_experiment-
metadata.yaml" -T 'gencode_mouse_m13' -C "$path/salmon_rna_seq.yaml" -I
"$path_to_fastq/$name-lec_1_r1.fq.gz" "$path_to_fastq/$name-
lec_1_r2.fq.gz" -O "$path/res/$name-lec_1" -S "$name-lec_1" -P '6' -X 'True'

```

The command is run by `Pypiper`: https://github.com/epigen/pypiper Inside of `salmon_rna_seq.py` I am also running commands to treat the data, the `salmon` part looks the following way:

```
def salmon(self):
    cmd = "salmon quant -i " + self.indexed_transcriptome
    cmd += " -l " + self.lib_type

    if self.trim:
        cmd += " -1 " + self.sample_files['r1']['fastq']['trimmed_paired']
        cmd += " -2 " + self.sample_files['r2']['fastq']['trimmed_paired']
    else:
        cmd += " -1 " + self.sample_files['r1']['fastq']['chastity']
        cmd += " -2 " + self.sample_files['r2']['fastq']['chastity']

    cmd += " -o " + self.sample_dirs['salmon']['base']
    cmd += " --numBootstraps=30"
self.pipe_manager.run(cmd, self.sample_files['salmon'])

```

Commands that the `Pypiper` runs are the following ones (`_commands.sh` `Pypiper` file):

```
zcat -c $PATH/res/old-lec_1/fastq/merged/old-lec_1_r1.fq.gz | grep -A 3
'^@.*[^:]*:N:[^:]*:' | grep -v '^\-\-$' | gzip -c > $PATH/res/old-lec_1/fastq
/chastity/old-lec_1_r1.fq.gz

zcat -c $PATH/res/old-lec_1/fastq/merged/old-lec_1_r2.fq.gz | grep -A 3
'^@.*[^:]*:N:[^:]*:' | grep -v '^\-\-$' | gzip -c > $PATH/res/old-lec_1/fastq
/chastity/old-lec_1_r2.fq.gz

java -jar $HOME/anaconda2/share/trimmomatic-0.36-5
/trimmomatic.jar PE -threads 6 -phred33 $PATH/res/old-lec_1/fastq/chastity/old-
lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r2.fq.gz $PATH/res/old-
lec_1/fastq/trimmed/old-lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/trimmed/unpaired
/old-lec_1_r1.fq.gz $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r2.fq.gz $PATH/res/old-
lec_1/fastq/trimmed/unpaired/old-lec_1_r2.fq.gz ILLUMINACLIP:$HOME/anaconda2/share/trimmomatic-0.36-5/adapters
/Small_List_Of_Adapters.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:15
MINLEN:36

fastqc $PATH/res/old-lec_1/fastq/merged/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/merged

fastqc $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/chastity

fastqc $PATH/res/old-lec_1/fastq/merged/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/merged

fastqc $PATH/res/old-lec_1/fastq/chastity/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/chastity

fastqc $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r1.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/trimmed

fastqc $PATH/res/old-lec_1/fastq/trimmed/old-lec_1_r2.fq.gz --outdir $PATH/res/old-
lec_1/fastqc/trimmed

salmon quant -i $HOME/Genomes/mouse/salmon
/gencode.vM13.transcripts.fa.sidx -l IU -1 $PATH/res/old-lec_1/fastq/trimmed
/old-lec_1_r1.fq.gz -2 $PATH/res/old-lec_1/fastq/trimmed/old-
lec_1_r2.fq.gz -o $PATH/res/old-lec_1/quantification --numBootstraps=30

```

`Edit`

After running manually `gzip -t` on a file:

```
-bash-4.2$ gzip -t old-lec_1_r1.fq.gz
gzip: old-lec_1_r1.fq.gz: unexpected end of file

```

So, I guess the file is corrupted and the issue is with some library that is generating the file. However, reinstalling manually other libraries do not help.

# Answer
The issue was caused by the corrupted input `.fastq` files that were damaged somehow upon uploading them from the local machine. We figured it out by `md5sum` command output comparison. We have not yet managed to actually run all of the script still, it is failing, but for another reason now.
> 1 votes

---
Tags: rna-seq, quantification, salmon
---


# Question
Title: Pipeline (Pypiper) returns error, but individual command works fine
I am running a script file on a cluster with `slurm` through `sbatch`. Locally it runs perfectly well on 2 different linux versions (`linux mint`, `ubuntu 12.04`), however, it returns `-9` error code while running on a cluster which is `CentOS 7.3`. When we are using `Pypiper` and run the whole `.py` script file, it fails when it starts running `trimmomatic` program

http://www.usadellab.org/cms/?page=trimmomatic

However, when running just `trimmomatic` program step by itself it does not fail and finishes successfully. We have a suspicion that something may be wrong with how `Pypiper` runs on the cluster, however, feeling rather stuck. Any suggestions would be greatly appreciated.

Here is the error that is printed out:

```
ILLUMINACLIP: Using 2 prefix pairs, 22 forward/reverse sequences, 0
forward only sequences, 0 reverse only sequences
</pre>
Process 109391 returned: (-9). Elapsed: 0:00:51. Peak memory: (Process:
5.211GB; Pipeline: 0.0GB)

### Pipeline failed at:  (12-14 15:50:47) elapsed:0:23:09 _TIME_

Total time: 0:23:09

Change status from running to failed
Traceback (most recent call last):
  File "...../salmon_rna_seq.py", line 344, in <module>
sys.exit(main())
File "...../salmon_rna_seq.py", line 304, in main
pipeline.process()
File "...../salmon_rna_seq.py", line 277, in process
self.trim_reads()
File "...../salmon_rna_seq.py", line 195, in trim_reads
self.pipe_manager.run(cmd, self.sample_files['r1']['fastq']['trimmed_paired'])
File "...../pypiper.py", line 508, in run
process_return_code, local_maxmem = self.callprint(cmd, shell, nofail, container)  # Run command
File "...../pypiper.py", line 691, in callprint
self.fail_pipeline(e)
File "...../pypiper.py", line 1111, in fail_pipeline
raise e
OSError: Subprocess returned nonzero result.

Pypiper terminating spawned child process 108475...
child process terminated
slurmstepd: error: Exceeded step memory limit at some point.

```

The last line may be informative. Although, we are still not sure what we could possibly do.

# Answer
Slurm is killing the process because it's using more memory per core than the default allocation (that's what `slurmstepd: error: Exceeded step memory limit at some point` means). I've never used pyPiper (like seemingly most people, I use snakemake), but presuming that you have it running either `srun` or `sbatch` at some point then specify something like `--mem-per-cpu 10G`.
> 3 votes

---
Tags: rna-seq
---


# Question
Title: changing color key range to specified range in heatmap.2 function
I have a tab separated text file as shown below

```
a   0.311   0.510   0.123   0.002   0.001   0.417   0.572   0.074   0.169
b   0.324   0.592   0.070   0.028   0.028   0.535   0.535   0.127   0.113

```

I am trying to use heatmap.2 function from gplots package and want to change the color key so that it ranges from 0 to 1

Following is my code:

```
d <- read.table("file.txt",sep="\t",header=F,row.names = 1)
colnames(d) <- c("-3","-2","-1","+1","+2","+3","+4","+5","+6")
d_matrix <- data.matrix(d)
colors<-seq(0,1,length.out=100)
heatmap.2(d_matrix,scale="none",col = redblue(75),trace="none",Colv = FALSE,Rowv = FALSE,breaks=colors)

```

But this gives me error:

```
Error in image.default(1:nc, 1:nr, x, xlim = 0.5 + c(0, nc), ylim = 0.5 +  :
  must have one more break than colour

```

I want to change to color key such that it ranges from 0 to 1. How can I do it in the function heatmap.2

I added breaks argument to the heatmap.2 function as suggested from this post here but no luck..

# Answer
You just need to add more colors. The following, for example, will work:

```
heatmap.2(d_matrix, scale="none", trace="none", Rowv=F, Colv=F, breaks=seq(0,1,0.01), dendrogram="none", col=redblue(100))

```

Note that `length(seq(0, 1, 0.01))` is 101, therefore `redblue()` must return one fewer (100) colors (or colours if you prefer extra vowels).
> 4 votes

---
Tags: r, bioconductor
---


# Question
Title: Convert population allele count into population allele frequency by using Perl or Python (translate available Bash into Perl or Python)
I already know the Bash AWK solution and also R code for the question that I am asking. But, my file is so big and with R it takes very long time and I am afraid of AWK's mistake for this computation, so I wonder if there is a Perl or Python solution for my problem.

I have a population allele count data, and it looks like this:

```
1   0   0   0   0   0   0   0   0   0   1   2   1   0   0   0   0
0   2   0   0   0   0   0   0   0   0   0   4   0   2   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0   0
2   2   0   0   2   1   0   0   0   0   2   4   0   0   0   2   0

```

Columns are populations and rows are SNPs. I have two rows per SNP (one row for the number of copies of allele "A" in each population and one row for the number of copies of allele "a").

In the example above, the first and second rows are number of alleles for the SNP1, the third and fourth rows are alleles for SNP 2, and like this until SNP20000000.

I want to calculate population allele frequency per SNP for all populations: frequency of allele A at SNP1 in population 1 = number of copies of allele A in population/total number of A/a gene copies in population and frequency of allele a at SNP1 in population 1 = number of copies of allele a in population/total number of A/a gene copies in the population. This means for each SNP number of copies of each allele should be divided to sum of the number of copies of allele "A" and "a" for each population. This is my desired output:

```
1   0   0   0   0   0   0   0   0   0   1   0.333333    1   0   0   0   0
0   1   0   0   0   0   0   0   0   0   0   0.666667    0   1   0   0   0
0   0   0   0   0   0   0   0   0   0   0   0           0   0   0   0   0
1   1   0   0   1   1   0   0   0   0   1   1           0   0   0   1   0

```

As I said, I have an R and Bash solution, but is there a way to do this estimation in Perl or Python?

Here is the Bash solution, but I do not know how to translate it into Perl or Python? I am learning Perl...

```
awk '{for(i = 1; i <= NF; i++)
    tally[i] += $i};
    (NR%2) == 1
    {for(i=1; i<=NF; i++)
        allele1[i] = $i};
        (NR%2) == 0
        {for(i = 1; i <= NF; i++)
        allele2[i] = $i;
            for(i = 1; i <= NF; i++)
                if(tally[i] == 0)
                    tally[i] = 1;
                for(i = 1; i <= NF; i++)
                    printf allele1[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    printf allele2[i]/tally[i]"\t";
                    printf "\n";
                for(i = 1; i <= NF; i++)
                    tally[i]=0}' MyData |
sed 's/\t$//g'

```

# Answer
I don't know how fast you need it to be, but this R solution runs on 20 million SNPs in under 2 min on my laptop:

```
allele_counts <-matrix(rep(c(1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 2, 1, 0, 0, 0, 0,
                     0, 2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 4, 0, 2, 0, 0, 0,
                     0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                     2, 2, 0, 0, 2, 1, 0, 0, 0, 0, 2, 4, 0, 0, 0, 2, 0), 10000000),
                   ncol=17, byrow=TRUE)
allele1<-allele_counts[ c(TRUE,FALSE), ]
allele2<-allele_counts[ c(FALSE,TRUE), ]
total_counts <- allele1+allele2
allele_counts[c(TRUE,FALSE), ] <-allele1/total_counts
allele_counts[c(FALSE,TRUE), ] <-allele2/total_counts

```

*Edit*

I decided to a head-to-head test of this solution vs. bli's pandas solution and your awk solution

I used readr for file IO:

```
library(readr)
allele_counts<-as.matrix(read_tsv("allele_counts.txt", col_names=F))
...
write_tsv(as.data.frame(allele_counts), "allele_freq_R.txt", col_names=F)

```

Here are the timings:

```
cat Benchmarks/R.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
103.7556        0:01:43 2211.16 2518.59 2209.11 2209.15 0.00    27.58   73.62

cat Benchmarks/Python.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
132.6425        0:02:12 4872.45 6042.28 4870.89 4870.92 2.67    46.59   99.92
cat Benchmarks/awk.txt
s       h:m:s   max_rss max_vms max_uss max_pss io_in   io_out  mean_load
126.4455        0:02:06 4.22    413.55  0.93    0.97    0.00    130.68  91.28

```

The R solution is slightly faster than the pandas solution and uses half as much memory. Your awk one-liner is about the same speed as the pandas script, but requires a tiny fraction of the memory footprint.
> 4 votes

# Answer
In python, if your files can fit in RAM, you might try using pandas and cytoolz.

```
import pandas as pd
from cytoolz import interleave

def odd(num):
    return num % 2

def even(num):
    return 1 - (num % 2)

# Read in two separate dataframes
alleles_1 = pd.read_table("snps.txt", skiprows=odd, sep="\s+", header=None)
alleles_2 = pd.read_table("snps.txt", skiprows=even, sep="\s+", header=None)
sums = alleles_1 + alleles_2
# https://stackoverflow.com/a/45565745/1878788
pd.DataFrame(list(interleave([
    alleles_1.div(sums).fillna(0).values,
    alleles_2.div(sums).fillna(0).values]))).to_csv(
        "freqs.txt", sep="\t", index=False, header=False)

```

I don't know how efficient this is. On my workstation it runs in 1m27 on data obtained by taking 1 million times the example given in the question. See heathobrien's answer for benchmarks on larger data.
> 0 votes

# Answer
Here is an implementation in pure python. Probably not the fastest, but you should not encounter any memory issues with this.

Note: this is python3, use as follows: `cat snp_counts.tsv | python allele_frequencies.py`

```
#!/usr/bin/env python
"""
allele_frequencies.py
"""
import sys

is_odd = True
for line in sys.stdin:
    parts = line.strip().split()
    if not parts:
        continue
    allele_counts = [ int(p) for p in parts ]
    if is_odd:
        allele1 = allele_counts
        is_odd = False
    else:
        allele2 = allele_counts
        allele_sum = [ sum([a1,a2]) for a1,a2 in zip(allele1,allele2) ]
        allele1_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele1, allele_sum) ]
        allele2_freq = [ a / tot if tot else 0.0 for a,tot in zip(allele2, allele_sum) ]
        print(*allele1_freq, sep = '\t')
        print(*allele2_freq, sep = '\t')
        is_odd = True

```
> 0 votes

# Answer
Here's my quick Perl attempt after misreading the question, which divides each row elements by the row total:

```
cat popCountData.txt | perl -lane '
  my $total=0;
  grep {$total += $_} @F; ## create per-row total
  foreach (@F) {printf("%0.2f ", ($total) ? $_ / $total : 0)} @F; ## divide elements by total
  print ""; ## print a line break (note: line will end with a space)
'
0.20 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.20 0.40 0.20 0.00 0.00 0.00 0.00
0.00 0.25 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.50 0.00 0.25 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.13 0.13 0.00 0.00 0.13 0.07 0.00 0.00 0.00 0.00 0.13 0.27 0.00 0.00 0.00 0.13 0.00

```

The challenge here is that the processing is in chunks of two lines, which means a couple of things:

1. The previous line (at least for even-numbered lines) needs to be stored
2. The sum of counts (or whatever statistic is necessary) will have different processing on each line

Unfortunately, I can't quite see how to get from your example input to your example output. For example, how does a count of 2 and 4 turn into a frequency of $1/3$ and $2/3$?
> -1 votes

---
Tags: python, perl
---


# Question
Title: Using Python, how to convert a pandas DataFrame into a VCF?
Let's say I have a pandas dataframe with fields `CHROM`, `POS`, `ALT`, `REF`. In this special case, I also wouldn't care about `ID`or `FILTER`, `INFO` could be blank (or meaningless) and we'll write `QUAL` as each 40.

```
import pandas as pd
example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)
print(df)

   CHROM      POS   REF ALT
0     20    14370     G   A
1     20    17330     T   A
2     20  1110696     A   G
3     20  1230237     T   C
4     20  1234567  GTCT   G

```

In this case, saving the output with the extension `.vcf` via `df.to_csv("myfile.vcf")` with the other required fields would created the tab-delimited format of the VCF, but it wouldn't include the mandatory header. https://samtools.github.io/hts-specs/VCFv4.1.pdf

```
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
...
#CHROM POS ID REF ALT QUAL FILTER INFO

```

The most popular package to work with VCF's appears to be `PyVCF`: https://pypi.python.org/pypi/PyVCF

However, one needs to begin with a VCF as a template, and add to this (from the documentation):

```
vcf_reader = vcf.Reader(filename='vcf/test/tb.vcf.gz')
vcf_writer = vcf.Writer(open('/dev/null', 'w'), vcf_reader)
for record in vcf_reader:
     vcf_writer.write_record(record)

```

It appears the only information necessary from `vcf.Reader(filename='vcf/test/tb.vcf.gz')` is the metadata?

```
print(vcf_reader.metadata)
OrderedDict([('fileformat', 'VCFv4.0'), ('fileDate', '20090805'), ('source', ['myImputationProgramV3.1']), ('reference', '1000GenomesPilot-NCBI36'), ('phasing', ['partial'])])

```

Otherwise, `tb.vcf.gz` already has variants, which I think means I would need to feed in a VCF with no variants?:

```
for record in vcf_reader:
     print(record)

Record(CHROM=20, POS=14370, REF=G, ALT=[A])
Record(CHROM=20, POS=17330, REF=T, ALT=[A])
Record(CHROM=20, POS=1110696, REF=A, ALT=[G, T])
Record(CHROM=20, POS=1230237, REF=T, ALT=[None])
Record(CHROM=20, POS=1234567, REF=GTCT, ALT=[G, GTACT])

```

# Answer
You will have to define the header either from an existing VCF file or hardcoded into you python script. Then write the header the the output VCF file then write the dataframe to the same file with the mode options set to 'a' to append to the end of the file.

You will need to order the dataframe columns to match the VCF spec and insert the empty/ default `ID`, `FILTER`, `INFO` and `QUAL` columns

```
import pandas as pd

example_dict = {'CHROM':[20, 20, 20, 20, 20], 'POS':[14370, 17330, 1110696, 1230237, 1234567], 'REF':['G', 'T', 'A', 'T', 'GTCT'], 'ALT':['A', 'A', 'G', 'C', 'G']}

df = pd.DataFrame(example_dict)

df['ID'] = ""
df['QUAL'] = 40
df['FILTER'] = ""
df['INFO'] = ""

df = df[['CHROM', 'POS', 'ID', 'REF', 'ALT', 'QUAL', 'FILTER', 'INFO']]

header = """##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
"""

output_VCF = "myfile.vcf"
with open(output_VCF, 'w') as vcf:
    vcf.write(header)

df.to_csv(output_VCF, sep="\t", mode='a', index=False)

```

Then `myfile.vcf` looks like

```
##fileformat=VCFv4.1
##fileDate=20090805
##source=myImputationProgramV3.1
##reference=file:///seq/references/
#CHROM POS ID REF ALT QUAL FILTER INFO
CHROM   POS ID  REF ALT QUAL    FILTER  INFO
20  14370       G   A   40
20  17330       T   A   40
20  1110696     A   G   40
20  1230237     T   C   40
20  1234567     GTCT    G   40

```
> 8 votes

---
Tags: python, vcf
---


# Question
Title: How is a principal component analysis for population structure made from a bam or other alignment file?
I'm curious how a PCA for population structure is built from an alignment/BAM or SAM file. How do programs like angsd and ngsCovar actually build a PCA? What exactly makes things come out as different and how can this be determined? If I have 3 alignments that look like this:

---

Locus1:

```
>Indiv1
ACGT
>Indiv2
ACGC
>Indiv3
ACGA

```

---

Locus2:

```
>Indiv1
GACT
>Indiv2
GACA
>Indiv3
GAGT

```

---

Locus3:

```
>Indiv1
GGTC
>Indiv2
GGCT
>Indiv3
GCTA

```

---

(These example alignments leave out quality scores for simplicity, assume that each has high quality across the board.)

What algorithm or system would be used to graph them on a PCA? How many dimensions would that PCA have? Is it possible to model this by hand so that I can understand what's going on/how it is built?

# Answer
> 4 votes
In general, the process works as such:

1. Variants are found at a number of loci in each sample (missing loci are often imputed).
2. Each sample is scored at each locus, often according to the number of copies of the reference allele that it has (so, 0-2).
3. The resulting matrix is used for PCA.

---
Tags: sam, sequence-alignment, algorithms, phylogenetics
---


# Question
Title: Finding a single open reading frame with ribosomal binding site, using Biopython
I'm given a Fasta file, containing a large DNA(over 115,000 long) sequence, and I am tasked with finding a **single** large open reading frame contained within the DNA sequence using Biopython.

I'm aware this has been asked before consistently however i can't find anything for actually finding the source of my ORF, but apologies for any overlap.

Now from other sources and the Biopython cookbook I've translated my sequence and found six open reading frames (three for each strand) and their positions within the sequence;

```
def find_orfs_with_trans(seq, trans_table, min_protein_length):
    answer = []
    seq_len = len(seq)
    for strand, nuc in [(+1, seq), (-1, seq.reverse_complement())]:
        for frame in range(3):
            trans = str(nuc[frame:].translate(trans_table))
            trans_len = len(trans)
            aa_start = 0
            aa_end = 0
            while aa_start < trans_len:
                aa_end = trans.find("*", aa_start)
                if aa_end == -1:
                    aa_end = trans_len
                if aa_end-aa_start >= min_protein_length:
                    if strand == 1:
                        start = frame+aa_start*3
                        end = min(seq_len,frame+aa_end*3+3)
                    else:
                        start = seq_len-frame-aa_end*3-3
                        end = seq_len-frame-aa_start*3
                    answer.append((start, end, strand,
                                   trans[aa_start:aa_end]))
                aa_start = aa_end+1
    answer.sort()
    return answer

orf_list = find_orfs_with_trans(record.seq, table, min_pro_len)

for start, end, strand, pro in orf_list:
print("%s...%s - length %i, strand %i, %i:%i" \
      % (pro[:30], pro[-3:], len(pro), strand, start, end))

```

**OUTPUT**

```
 TNRQVYGGTLQSLRTGTGIYSRLASSPTNR...CEA - length 160, strand 1, 7950:8433
 GIPPGRTEGLGRYVHGESIFLEATLVPEPQ...VQS - length 171, strand -1, 19275:19791
 ISVHLRRYFSVLLRAPVALNADSRVAGPLD...ITP - length 190, strand 1, 34079:34652
 GCVNGNFPDYRVAVDDPGALVVGGELGQAL...VNT - length 771, strand -1, 39335:41651
 LLKLLASRLTWLLVLSAALGFLTSVSYRLG...SPG - length 235, strand -1, 39358:40066
 PGDVIVSKVPADNLRPRMSEINDFTNSIII...VLV - length 826, strand 1, 39362:41843
 GPGRQSSSAYERDQRLHQQHHHQRERHQVR...QRS - length 764, strand 1, 39385:41680
 LLLSSTSSASCWLKKFFFRLLLLSRSVALL...LAL - length 158, strand -1, 40459:40936
 TRLIIRPPSHTTLYYDPYPQPLKVSSSVPC...RIT - length 159, strand -1, 54663:55143
 IGAFSRGYNQGWVLIEAQEVRTTSPLFPRS...RDQ - length 201, strand 1, 56104:56710
 KARGSSLGAIYCRRLKRRSLPSWQIESFNC...TNS - length 150, strand -1, 62497:62950
 RRPARAGAELNKDDIRDTCLLSCSEVDNKV...VTL - length 168, strand 1, 72339:72846
 SGPTTVRTSAPQRCWRTSLKPKLSLSDSTE...PVA - length 177, strand -1, 114494:115028
 ERELWLETGPPTPLWGTGPDCSRAALNRQR...PRR - length 157, strand 1, 114953:115427

```

However I'm unsure if these open reading frames contain the ribosome binding site(consensus sequence of *AAGGAGGTG* and occurring between **6 to 9** nucleotides upstream from the positive 5' strand). Is there anyway to add to my code to include this? or would it already be included with Biopython?

My specification says to find **the** single open reading frame however since none of these begin the the *M* start codon i'm a bit lost how they can be so, I have no information on where my sequence is from (something I have to find out) so using the standard **translation table 1**.

# Answer
You should start by finding all instances of `AAGGAGGTG` and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.

There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.
> 2 votes

# Answer
From what I can tell, the program output does not indicate whether the open reading frames (ORFs) are preceded by ribosome binding sites (RBSs). I'm not sure what parameters are used when establishing ORFs. To verify your RBS requirement, index the input DNA (using given ORF locations from biopython output, pay attention to strands, which are designated by 1/-1) and search upstream for a sequence similar to your consensus. You'd have to pre-define your criteria for similarity to consensus in order to find it programmatically (or you can do it visually to better avoid false negatives). Perhaps try a search that permits ambiguous/non-matching third nucleotides in codons (like a regular expression such as "AA.GA.GT." where "." can be any nucleotide). Biologically, third nucleotide changes are least probable to drastically alter the translated residue (more likely to encode amino acid with similar chemical properties).
> 2 votes

---
Tags: python, homework
---


title: Normalizing RNAseq for PCA and CCA

text:
# Question
Title: Normalizing RNAseq for PCA and CCA
Usually the expression data is transformed to log space using either RPKM, FPKM or CPM, this is required when looking for differential expression because the data is tested against the normal distribution(limma) or the negative bionimal distribution (DESeq2).

In weighted gene co-expression analysis (WGCNA) there is also the recommendation to normalize the data (see FAQ 4) but says "RPKM, FPKM, or simply normalized counts doesn't make a whole lot of difference for WGCNA" because it is based in correlations.

For PCAs the data should be on the same scale to avoid missing the scale as a factor. See this for more references. But should it be normalized too to look like a normal distribution using logarithms?

How does this extend to other methods like CCA? Should the data look like a normal distribution for CCA? As it is based on correlation I don't expect much changes, so I am using raw counts (RSEM output), but the examples I've seen use normalized data.

This question comes after reading a comment

# Answer
I've given an answer about using how I normalise for PCA here. PCA works best if the data are normal-ish, but as models go it is fairly robust to fairly large deviations from normal. Transforming is not essential, but if data can be easily transformed into something that looks a bit more normal, then conclusions derived from looking at the PCA should be a bit more robust.

I'm not familiar with canonical-correlation analysis (CCA), but because I see sums and square roots in the formulas, I would assume it has similar data cleaning requirements to PCA. Unfortunately the R documentation for `cancor` doesn't have any information about any assumptions made on the data.
> 2 votes

---
Tags: rna-seq, normalization
---


# Question
Title: How to compute the chance of failing to detect a gene given the detection limit of a protocol
In Shapiro et al., when discussing about loss of molecules as source of error in single-cell sequencing, it is written that:

> Another source of error is losses, which can be severe. The detection limit of published protocols is $5$–$10$ molecules of mRNA. If, as seems likely, the limit of detection is primarily determined by losses during sample preparation, this would indicate that $80$–$90\\%$ of mRNA was lost. Or, to put it the other way around, a $90\\%$ loss leads to an approximately $50\\%$ chance of failing to detect a gene that is expressed at a level of seven mRNA molecules (from the binomial distribution).

How is this probability computed using the binomial distribution? I thought that $90\\%$ loss corresponds to $5$ detected molecules, and I assume that $k=7$ for the binomial calculation, but I am unable to go further.

# Answer
> 3 votes
A 90% loss can be rephrased as a 10% chance of detecting anything. So what we want to find is the probability of detecting 0 molecules, when we start with 7 and have 10% probability of success. Once can do that in R as follows:

```
> pbinom(0, 7, 0.1)
0.4782969

```

So ~50%, as they stated. I suspect that part of the confusion arises from the fact that the detection limit is due to lowly-expressed genes/transcript being heavily affected by this loss. So the probability of detecting a single molecule out of 4 original molecules (assuming 90% mRNA loss) is 34%, for 3 molecules it's 27%, for 2 it's 19% and for 1 it's 10%. I think the threshold of 5 is mentioned more because it's a nice round number than there's anything particularly different in detecting a gene with 4 vs. 5 molecules (34 vs. 41% probability).

---
Tags: rna-seq, rna, single-cell
---


# Question
Title: Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions
I am having a `SingleCellExperiment` (sce) object:

> sce \<- SingleCellExperiment(list(counts=UMI\_count))

which has two dimensions:

```
dim(sce)
[1] 46761 68027

```

When I am trying to get column sums for all rows the following does not work:

```
colSums(assay(sce), na.rm = FALSE, dims = 1)

```

Giving me an error:

> Error in colSums(assay(sce)) : 'x' must be an array of at least two dimensions

Apparently the issue was already described:

https://support.bioconductor.org/p/100954/

However, I can not find the solution. Is there another way of getting sums of columns of either `SingleCellExperiment` object or just `dgcMatrix`?

# Answer
> 0 votes
Ultimately I just implemented it manually with the precursor to `SingleCellExperiment` matrix (`UMI_count` which is `dgcMatrix` type) object in the following way:

```
# Removing genes with 0 - counts in all of the columns
num_removed = 0
for(i in 1:nrow(UMI_count)) {
  idx = i - num_removed
  if(sum(UMI_count[idx,]) < 1) {
    UMI_count <- UMI_count[-idx,]
    num_removed = num_removed + 1
    print(i)
  }
}

sce <- SingleCellExperiment(list(counts=UMI_count))

```

I am getting sums of the columns for each row and then filtering out rows with the `sum < 1`. So, `sum` function for the slices - rows (`numeric` type) is working.

---
Tags: rna-seq, r
---


# Question
Title: calculating nucleotide frequency per column
I have some sequences shown below

```
CAGGTAGCC
CCGGTCAGA
AGGGTTTGA
TTGGTGAGG
CAAGTATGA
ACTGTATGC
CTGGTAACC
TATGTACTG
GCTGTGAGA
CAGGTGGGC
TCAGTGAGA
GGGGTGAGT
TGGGTATGT
GAGGTGAGA
CAGGTGGAG

```

Each line has 9 nucleotides. Consider it to be 9 columns.I want to calculate the nucleotide frequency of each nucleotide for each of the 9 columns. For example 1st column will have these bases C,C,A,T,C,A etc Out put should be something like this

```
A   0.25    0.34    0.56    0.43    0.00    0.90    0.45    0.34    0.31
C   0.45    0.40    0.90    0.00    0.40    0.90    0.30    0.25    0.2
G   0.00    0.00    0.34    1.00    0.30    0.30    0.35    0.90    0.1
T   0.24    0.56    0.00    0.00    1.00    0.34    0.45    0.35    0.36

```

---

**Note**, I just made up the numbers to show you the output file format

# Answer
For anyone interested in doing this from any sort of alignment file, I've implemented a position frequency matrix function in AlignBuddy.

input:

```
 15 9
seq_1      CAGGTAGCC
seq_2      CCGGTCAGA
seq_3      AGGGTTTGA
seq_4      TTGGTGAGG
seq_5      CAAGTATGA
seq_6      ACTGTATGC
seq_7      CTGGTAACC
seq_8      TATGTACTG
seq_9      GCTGTGAGA
seq_10     CAGGTGGGC
seq_11     TCAGTGAGA
seq_12     GGGGTGAGT
seq_13     TGGGTATGT
seq_14     GAGGTGAGA
seq_15     CAGGTGGAG

```

Command:

```
:$ alignbuddy input.phy --pos_freq_mat

```

Output:

```
### Alignment 1 ###
A       0.133   0.400   0.133   0.000   0.000   0.400   0.467   0.067   0.400
C       0.400   0.267   0.000   0.000   0.000   0.067   0.067   0.133   0.267
G       0.200   0.200   0.667   1.000   0.000   0.467   0.200   0.733   0.200
T       0.267   0.133   0.200   0.000   1.000   0.067   0.267   0.067   0.133

```
> 3 votes

# Answer
awk

```
awk '{L=length($1);for(i=1;i<=L;i++) {B=substr($1,i,1);T[i][B]++;}} END{for(BI=0;BI<4;BI++) {B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));printf("%s",B); for(i in T) {tot=0.0;for(B2 in T[i]){tot+=T[i][B2];}printf("\t%0.2f",(T[i][B]/tot));} printf("\n");}}' input.txt

A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

```

Or, expanded for clarity:

```
awk '{
        L=length($1);
        for(i=1;i<=L;i++) {
            B=substr($1,i,1);
            T[i][B]++;
        }
      }
      END{
            for(BI=0;BI<4;BI++) {
                B=(BI==0?"A":(BI==1?"C":(BI==2?"G":"T")));
                printf("%s",B);
                for(i in T) {
                    tot=0.0;
                    for(B2 in T[i]){
                        tot+=T[i][B2];
                    }
                printf("\t%0.2f",(T[i][B]/tot));
                }
            printf("\n");
            }
}' input.txt

```
> 5 votes

# Answer
Here is one example of how to do this with a bit of python. Alternatively one could create strings of each column and using `letterFrequency()` from the `Biostrings` package.

```
#Make a list of hashes
hl = []
for i in range(9):
    hl.append({'A': 0, 'C': 0, 'G': 0, 'T': 0})

f = open("foo.txt")  # CHANGE ME
nLines = 0
for line in f:
    for idx, c in enumerate(line.strip()):
        hl[idx][c] += 1
    nLines += 1
f.close()

nLines = float(nLines)
for char in ['A', 'C', 'G', 'T']:
    print("{}\t{}".format(char, "\t".join(["{:0.2f}".format(x[char]/nLines) for x in hl])))

```

The output of your example is then:

```
A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

```
> 4 votes

# Answer
Here's a Perl approach:

```
#!/usr/bin/env perl
use strict;
my %counts;

## Read the input file line by line
while (my $line = <>) {
  print;
  ## remove trailing '\n' characters
  chomp $line;
  ## split the line into an array at every character
  my @columns=split(/./, $line);
  ## iterate over the array from the first to the last position
  for my $i (0..$#columns){
    ## The nucleotide at this position
    my $nt = $columns[$i];
    ## Save the count in this hash of arrays. The keys are the
    ## nucleotides, and the value at each position is increased
    ## when that nucleotide is found in that position.
    $counts{$nt}[$i]++;
  }
}
## Iterate over the counts hash
for my $nt (sort keys(%counts)){
  print "$nt\t";
  ## dereference the array stored in the hash
  my @countsForThisNt = @{$counts{$nt}};
  ## iterate over the counts for each position for this nt
  for (my $l=0;$l<=$#countsForThisNt; $l++) {#
    ## If the value for this position isn't defined,
    ## set it to 0.
    $countsForThisNt[$l]||=0;
    ## Print all the things
    printf "%.2f\t", $countsForThisNt[$l]/$., $l;
  }
  print "\n";
}

```

Save the script somewhere in your PATH, make it executable and run:

```
$ foo.pl file
A   0.13    0.40    0.13    0.00    0.00    0.40    0.47    0.07    0.40
C   0.40    0.27    0.00    0.00    0.00    0.07    0.07    0.13    0.27
G   0.20    0.20    0.67    1.00    0.00    0.47    0.20    0.73    0.20
T   0.27    0.13    0.20    0.00    1.00    0.07    0.27    0.07    0.13

```

Alternatively, if you're into the whole brevity thing, and enjoy some golfing, here's the same thing as a one-liner:

```
 perl -ne 'chomp;@F=split(//);for$i(0..$#F){$k{$F[$i]}[$i]++}}{for $nt(sort keys(%k)){print"$nt\t";for$i(0..$#{$k{$nt}}){$g=$k{$nt}[$i]||0; printf"%.2f\t",$g/$.}print"\n";}' file

```
> 3 votes

# Answer
Here an example in R using `Biostrings` and `letterFrequency` (as suggested by Devon Ryan).

```
library(Biostrings)

data <- read.table("DNA.txt", stringsAsFactors = F)

new <- matrix(nrow = 9, ncol = 15)

for(i in 1:9){
  for(j in 1:15){
    new[i,j] <- substring(data[j,], i, i)
  }
}

countTable <- matrix(nrow = 9, ncol = 4)
for(i in 1:9){
  columnSeq <- DNAStringSet(paste0(new[i,], collapse = ""))
  columnCounts <- letterFrequency(columnSeq, letters = "ACGT", OR = 0)
  countTable[i,] <- columnCounts
}
colnames(countTable) <- c("A", "C", "G", "T")

freqTable <- countTable/15
round(t(freqTable), digit = 2)
  [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9]
A 0.13 0.40 0.13    0    0 0.40 0.47 0.07 0.40
C 0.40 0.27 0.00    0    0 0.07 0.07 0.13 0.27
G 0.20 0.20 0.67    1    0 0.47 0.20 0.73 0.20
T 0.27 0.13 0.20    0    1 0.07 0.27 0.07 0.13

```
> 2 votes

# Answer
benn’s answer works but is very un-R-like: you wouldn’t normally iteratively assign to items in a matrix: R has powerful vectorised operations which make this manual work unnecessary. Here’s a more idiomatic solution:

```
sequences = readLines('sequences.txt')
bases_matrix = do.call(rbind, strsplit(sequences, ''))

apply(bases_matrix, 2L, function (col) {
    str = DNAString(paste(col, collapse = ''))
    letterFrequency(str, letters = 'ACGT', OR = 0L, as.prob = TRUE)
})

```

This uses the same Bioconductor packages. Since this is such a simple problem, it can also be written without Bioconductor:

```
bases = strsplit(sequences, '')
# Use a data.frame here so we can use factors in the next step:
# R does not support matrices of factors. Ugh.
bases_by_column = setNames(do.call(rbind.data.frame, bases), seq_along(bases[[1L]]))
# Ensure that every column will be a complete set of ACGT frequencies
bases_by_column = lapply(bases_by_column, factor, c('A', 'C', 'G', 'T'))
sapply(lapply(bases_by_columns, table), prop.table)

```

Using modern R idioms from the ‘magrittr’ package, I’d write this as a pipeline; this very directly shows the sequence of transformations.

```
do.call(rbind.data.frame, bases) %>%
    setNames(seq_along(bases[[1L]])) %>%
    lapply(factor, c('A', 'C', 'G', 'T')) %>%
    lapply(table) %>%
    sapply(prop.table)

```
> 2 votes

# Answer
Here is how you get the distributions by column using the `TraMineR` R package.

```
library(TraMineR)

sts.data <- c(
  "CAGGTAGCC",
  "CCGGTCAGA",
  "AGGGTTTGA",
  "TTGGTGAGG",
  "CAAGTATGA",
  "ACTGTATGC",
  "CTGGTAACC",
  "TATGTACTG",
  "GCTGTGAGA",
  "CAGGTGGGC",
  "TCAGTGAGA",
  "GGGGTGAGT",
  "TGGGTATGT",
  "GAGGTGAGA",
  "CAGGTGGAG"
)
seq <- seqdef(seqdecomp(sts.data,sep=''))

seqstatd(seq)

```

and here is the outcome of the `seqstatd` function

```
  [State frequencies]
   [1]  [2]  [3] [4] [5]   [6]   [7]   [8]  [9]
A 0.13 0.40 0.13   0   0 0.400 0.467 0.067 0.40
C 0.40 0.27 0.00   0   0 0.067 0.067 0.133 0.27
G 0.20 0.20 0.67   1   0 0.467 0.200 0.733 0.20
T 0.27 0.13 0.20   0   1 0.067 0.267 0.067 0.13

   [Valid states]
  [1] [2] [3] [4] [5] [6] [7] [8] [9]
N  15  15  15  15  15  15  15  15  15

   [Entropy index]
   [1]  [2]  [3] [4] [5]  [6]  [7]  [8]  [9]
H 0.94 0.94 0.62   0   0 0.78 0.87 0.62 0.94

```
> 1 votes

---
Tags: sequence-analysis, nucleotide-models
---


# Question
Title: What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?
What software/tool can I use to visualise a multiple sequence alignment in the xmfa format as produced by command-line tool Mauve?

I am particularly interested in seeing the single nucleotide polymorphisms on the screen.

I am doing a multiple sequence alignment of about a dozen whole genomes, all in a single contig.

Link to Mauve

# Answer
> 2 votes
Typically I use the Mauve viewer to visualize the alignments. Are there things you find lacking in that viewer?

There are multiple threads about converting xmfa to other formats

# Answer
> 1 votes
I supervised a student a while ago, and she made a tool to focus on (i.e., visualize) SNPs in multiple alignments. The tool is called ADOMA, and is quite simple. It uses clustal for alignment, so I am not sure it if it may work for your Mauve data. If your data are xmfa, you'll probably have to make fasta of it first.

So I am not sure if this is a solution, but maybe the tool can be of use for your research.

---
Tags: sequence-alignment, software-recommendation
---


# Question
Title: Getting structure and sequence related to PDB IDs
I have two kinds of interactions: transient and stable. We are supposed to work on stable interaction, like interactions between two monomers in a heterodimer.

In a heterodimer there are two chains in which there are some residues between monomers which take part in the interaction of monomers and building the structure of heterodimer, and there are other residues interacting inside each monomer for building the structure of that monomer. There are some PDB IDs related to heterodimers in PDB.

Now imagine that we are trying to show which residues are physically interacting between two monomers in a heterodimer, and then we are trying two show the site of these physically interacting residues in the sequence related to that special structure. For that, we need to find the sequence related to that structure and we can download it from PDB.

But if we want to make a multiple sequence alignment of that sequence, first we need to find the refseq sequence of that and then we must blast the sequence and find homologues of that sequence.

Here there are some challenges:

1. we must find residues interacting in the structure of heterodimer between two monomers.
2. we must find the sequences of monomer chains of the heterodimer and map the interacting residues in the structure to the sequence.
3. we know that when a structure is solved and its monomers are sequenced, may be they are not able to completely sequence the heterodimer and residue numbers are maybe different from the related refseq sequence.

I would like to know how can we find the related refseq sequence (or sequences, in the case of heterodimer) and then how can we map the physically interacting residues in the structure of PDB ID to refseq sequence (finding the sites of these interacting residues in that refseq sequence) regarding three challenges mentioned above.

Imagine that I have 100 PDB IDs of 100 nonredundant heterodimers and I would like to find and download the structures, sequence and finally refseq sequence. What should I do?

# Answer
> 9 votes
You can download seqs and structures based on a list of PDB ids using http://www.rcsb.org/pdb/download/download.do#FASTA

# Answer
> 2 votes
The post might seem old, but I wish to share my thaughts based on how I solved similar problem. Here it goes:

I will use an example to illustrate. Other posts have already explained how you could get the protein ID from PDB and so on. I used the database PDBSUM (www.ebi.ac.uk/pdbsum/) and with protein PDB ID: 1KJF.

By following these steps I got to the results below: Go to --\> www.ebi.ac.uk/pdbsum/ --\>(put in PDB ID) or (others e.g UniProt ID) --\> Enter --\> Click on "Ligands" --\>(a LIGPLOT appears at the botoom of the page). You can as well search just by pasting in the sequence.

Results: You can see that 1KJF is made up of 3 chains A, B & P and both chains have clearly anntated and interacting residues. You can now see clearly read inidividual types of interactions.

Another good option is http://protein.bio.unipd.it/ring. It shows you a graph plot of all inteacting residues and how they are linked to each other. You can adjust the default result settings to meet your needs.

Hope it helps. :)

---
Tags: protein-structure
---


# Question
Title: R: 'Matrix' can not be unloaded, but 'writeMM' method not found
I want to save a sparse and very large `dgcMatrix` onto disc. I read that it could be done with `writeMM` method. So, when I am trying to do that:

```
writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset")

```

`Rstudio` is giving me an error:

```
Error in writeMM(UMI_count, "/home/gene_count_filtered/filtered_dataset") :
could not find function "writeMM"

```

Ok, I read that this function is in the `Matrix` package. I installed it, but when I try loading it:

```
library(Matrix)

```

I am getting another error due to other libraries using `Matrix` already on the lower level, as I understand:

```
Error in value[[3L]](cond) :
Package ‘Matrix’ version 1.2.11 cannot be unloaded:
Error in unloadNamespace(package) : namespace ‘Matrix’ is imported by
‘scran’, ‘scater’, ‘SummarizedExperiment’ so cannot be unloaded

```

I also tried just`write.matrix(UMI_count, "/home/nikita/filtered_dataset")`, but got an error:

> Cholmod error 'out of memory' at file ../Core/cholmod\_memory.c, line 147

What could be the solution to the problem? I just need to somehow save a very large `dgcMatrix` on disc.

# Answer
Ultimately I just restarted `Rstudio` as `Konrad Rudolph` suggested, imported `Matrix` library, and then `writeMM` method worked just fine. However, `rownames` were converted into just numbers, but I needed to keep the gene names instead as it was in `dgcMatrix`, so after some search I found a different way of doing that from here:

http://www.sthda.com/english/wiki/saving-data-into-r-data-format-rds-and-rdata

So, I just did `save()` (and then, to load back - `load()`) and the `rownames` are maintained like I wanted it to be. So, storing it now as `.RData` file.
> 1 votes

---
Tags: rna-seq, r
---


# Question
Title: How to remove all BAM read groups from all reads (not just the header)?
I have problem with one my BAMs---it appears to have invalid read groups.

Normally when I have such a problem, I remove all the read groups from the BAM header as follows:

```
samtools view -H your.bam | grep -v "^@RG" | samtools reheader - original.bam > new.bam

```

and then I would execute `AddOrReplaceReadGroups` from Picard:

http://broadinstitute.github.io/picard/command-line-overview.html#AddOrReplaceReadGroups

I have tried this, and there must still be reads with a read group in the optional tags `RG:Z:`.

How to do I definitively remove all read groups from a BAM, including all of the reads?

# Answer
You might have to manually strip those auxiliary tags off:

```
samtools view -h your.bam | grep -v "^@RG" | sed "s/\tRG:Z:[^\t]*//" | samtools view -bo your_fixed.bam -

```

The `sed` bit is searching for the aux tag and removing everything up to the next tab.
> 7 votes

---
Tags: sam, samtools, reads, picard
---


# Question
Title: Using column 2 of one file to match with two columns of another file, and append
I have file 1 like following:

---

```
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

```

And file 2 as following

---

```
1   15776231    2   5008    G:5002  A:6 1   16855677    2   5008    A:5003  C:5

```

I am trying to write a code(but failing!) to do the following:

if `$2` of file2 comes in between `$2` and `$3` of file1 (between the interval of `$2` and `$3`) and if `$8` of file2 comes in between `$5` and `$6` of file1 (between the interval of `$5` and `$6`), then the output would be `$1,$2,$3,$4` from file1 and `$2,$3,$4,$5,$6` from file2 and `$5,$6,$7` from file1 and `$8,$9,$10,$11,$12` from file 2, all in one line. So something like this:

```
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

```

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2. Is there any way to make bedtools work on this?

(There could be many lines from file2 intersecting between the same line in file1,so that line could be repeated)

# Answer
Pipe a modified form of the second file to BEDOPS `bedmap` and the first file, then pipe *that* result to cut out desired columns:

```
$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

```

This should run pretty quickly and use very little memory, as it uses Unix streams on (presumably) sorted BED data.

You'll want to test the column indices in my answer by running this on your test inputs, but that's an easy adjustment for the larger dataset, in case I missed a column.

If you're not sure about sort order, pipe to `sort-bed`:

```
$ sort-bed first.unsorted.bed > first.bed
$ awk -v OFS='\t' '{ print $1, $2, ($2+1), $3, $4, $5, $6, $7, $8, $9, $10; }' second.txt | sort-bed - | bedmap --echo --echo-map first.bed - | cut -f1-4,8,9,11- - > answer.bed

```

Use streams where you can. It's a different way of thinking but can pay huge dividends.
> 3 votes

# Answer
Here's one approach in Perl. Note that I have not tested this extensively, and I am not confident it works. Please test first and let me know if there are problems.

```
#!/usr/bin/env perl
use strict;
use warnings;
use 5.10.0;
my (@file1,%file2);
open(my $f1, '-|',"sort -nk2 $ARGV[0]");
while (<$f1>) {
  chomp;
  my @fields = split(/\s+/);
  push @file1, [ @fields ];
}
close($f1);

open(my $f2, "sort -nk2 $ARGV[1] |");
line:while (<$f2>) {
  chomp;
  my @fields = split(/\s+/);
  foreach my $file1Line (@file1) {
    my $file1Start1 = $file1Line->[1];
    my $file1End1 = $file1Line->[2];
    my $file1Start2 = $file1Line->[4];
    my $file1End2 = $file1Line->[5];
    my $file2Start1 = $fields[1];
    my $file2End1 = $fields[2];
    my $file2Start2 = $fields[7];
    my $file2End2 = $fields[8];
    my @file1Fields = @$file1Line;
    if ($file1Start1 < $file2Start1 &&
        $file1End1 > $file2End1 &&
        $file1Start2 < $file2Start2 &&
        $file1End2 > $file2End2) {
      say join "\t", @{$file1Line}[0..3], @fields[1..5], @{$file1Line}[4..6], @fields[7..11];
    }
    elsif ($file1Start1 > $file2Start1) {
      last;
    }
  }
}
close($f2);

```

Save the script as `foo.pl` somewhere in your `$PATH`, make it executable and run as:

```
$ foo.pl file1.bed file2.bed
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776231    2   5008    G:5002  A:6 16855676    16855696    CTGTCCAGCAGAGGGCGGTG    16855677    2   5008    A:5003  C:5

```
> 1 votes

---
Tags: vcf, bedtools, bioawk, text-processing
---


# Question
Title: replacing SNPs with missing calls with a specific string
I have a big file containing 27 columns and nearly 6 million rows. The following is a little example of my file

```
head data
0.65   0.722222   1.0      0.75     0
0.35   0.277778   0.0      0.25     0
0      0.666667   0.75     0.5      0.5625
0      0.333333   0.25     0.5      0.4375

```

Rows are SNPs and columns are individuals. and I have "2 rows per each SNP" (one for reference allele and the other for alternative allele). In the example above I showed data for 2 SNPs (rows 1 and 2 correspond to SNP1 and rows 3 and 4 correspond to SNP2). I want to check if for each SNP, both ref and alt allele frequencies are 0, replace them with 9. this is my desired output:

```
head desired
0.65   0.722222   1.0      0.75     9
0.35   0.277778   0.0      0.25     9
9      0.666667   0.75     0.5      0.5625
9      0.333333   0.25     0.5      0.4375

```

# Answer
A possible approach in plain python:

```
#!/usr/bin/env python3

import sys

def both_zero(freq1, freq2):
    return freq1 == "0" and freq2 == "0"

with open(sys.argv[1], "r") as snp_file:
    for line1 in snp_file:
        # Consume one extra line to get the frequencies
        # for the alternative allele
        line2 = snp_file.readline()
        # "zip" the frequencies in pairs
        freq_pairs = zip(line1.strip().split(), line2.strip().split())
        recoded_pairs = [
            ("9", "9") if both_zero(freq1, freq2) else (freq1, freq2)
            for (freq1, freq2) in freq_pairs]
        # "unzip" the recoded frequencies into separate objects
        # (see https://stackoverflow.com/a/19343/1878788)
        (recoded1, recoded2) = zip(*recoded_pairs)
        print(*recoded1, sep="\t")
        print(*recoded2, sep="\t")

```

Calling it with your example snp frequencies file as command-line argument:

```
$ ./recode_absent.py snps.txt
0.65    0.722222        1.0     0.75    9
0.35    0.277778        0.0     0.25    9
9       0.666667        0.75    0.5     0.5625
9       0.333333        0.25    0.5     0.4375

```
> 2 votes

---
Tags: python, perl, text-processing
---


# Question
Title: I am getting this error when running a test data set for Lyve-SET on RHEL 7. Is this an MPileup error or something else?
This is the error that I am getting:

```
launch_varscan.pl: Creating a pileup /home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.mpileup
removed ‘/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.1239.mpileup’
launch_varscan.pl: Running mpileup:
  echo "NC_001416.1
" | xargs -P 1 -n 1 -I {} sh -c 'echo "MPileup on {}" >&2; samtools mpileup -f lambda/reference/reference.fasta
-Q 0 -B  --region "{}" /home/user/lyve-SET-1.1.4f/lambda/bam/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam >
/home/user/lyve-SET-1.1.4f/lambda/tmp/NC001416.fasta.wgsim.fastq.gz-reference.sorted.bam.$$.mpileup'
MPileup on NC_001416.1
mpileup: unrecognized option '--region'
open: No such file or directory
[mpileup] failed to open NC_001416.1: No such file or directory
launch_varscan.pl: main::mpileup: ERROR with xargs and samtools mpileup
Stopped at /home/arahat/lyve-SET-1.1.4f/scripts/launch_varscan.pl line 106.
QSUB ERROR
31488
launch_set.pl: Schedule::SGELK::command: ERROR with command: Inappropriate ioctl for device
  /usr/bin/perl /home/user/lyve-SET-1.1.4f/lambda/log/SGELK/qsub.2705.pl
Stopped at /home/user/lyve-SET-1.1.4f/scripts/../lib/Schedule/SGELK.pm line 661.
set_test.pl: ERROR: died while trying to run COMMAND
  launch_set.pl lambda --numcpus 1 --numnodes 20
set_test.pl: See if you have the correct software installed for this command.
6400 at /home/arahat/lyve-SET-1.1.4f/scripts/set_test.pl line 94.
make: *** [test] Error 25

```

I think the main problem is the "mpileup: unrecognized option --region", which seems to be blank. In the beginning of the block, it says creating pileup and then it immediately removes the file, probably to create the pipeline file.

I cannot change the code for Lyve-SET pipelining or create shell scripts for the pipeline because I am only configuring my server for use with the genome sequencing pipeline. Would this be a problem with the data set or some tool that I need to install? Any insight would be appreciated.

# Answer
The samtools v0.1.19 manual confirms that there is no `--region` option. In that version there is only the `-r` argument to pass in the targeted region. If you update to a newer version of samtools this should execute correctly
> 0 votes

---
Tags: samtools, mpileup
---


# Question
Title: Which identifier do I have to use when I want to add a ##Fasta section to a gff?
I have a gff file like this

```
scaffold_1  source  mRNA    2987526 2992430 .   -   .   ID=protein_id_68892;Name="foobar";locus_tag=my_organism_68892;translation=MHTGDALEGSTGNVSILV...
scaffold_1  source  exon    2987526 2987805 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2987526 2987805 .   -   1   name "foobar"; proteinId 68892; exonNumber 5;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2988610 2988745 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2988610 2988745 .   -   2   name "foobar"; proteinId 68892; exonNumber 4;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2989157 2989645 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2989157 2989645 .   -   2   name "foobar"; proteinId 68892; exonNumber 3;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2990903 2991541 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2990903 2991541 .   -   2   name "foobar"; proteinId 68892; exonNumber 2;locus_tag=my_organism_68892;Parent=protein_id_68892
scaffold_1  source  exon    2991587 2992430 .   -   .   name "foobar"; transcriptId 68892;locus_tag=my_organism_68892
scaffold_1  source  CDS     2991587 2992430 .   -   0   name "foobar"; proteinId 68892; exonNumber 1;locus_tag=my_organism_68892;Parent=protein_id_68892

```

The `mRNA`entry was not present at first and I added it in. I also have the protein fast file available and would now like to add the amino acid sequences to the gff. I originally added them in column 9 to the `mRNA`entry which is apparently incorrect. Here it is written that the sequences should go at the end of the file in the `##Fasta`section. However, I struggle to figure out which identifier I should use to link it to my `mRNA` entries.

Do I have to use

```
##Fasta
>protein_id_68892
MHTGDALEGSTGNVSILV...

```

or does it work with the locus\_tag

```
##Fasta
>my_organism_68892
MHTGDALEGSTGNVSILV...

```

or is it something completely else? Where and how do can I add the translation information?

In the link above (and also in the comments it is suggested) the identifier in the first column is used, in my case that would be `scaffold_1`. However, I do have a lot of entries with `scaffold_1`in my file so a link between `scaffold_1`and an amino acid sequence would not be unique as I would have

```
##Fasta
>scaffold_1
MHTGDALEGSTGNVSILV...
>scaffold_1
GDAHTEGSTGNVSILV...

```

. In the end, I want to create a genbank file (which does not exist yet) where I can link a CDS to an amino acid sequence.

# Answer
> 2 votes
Unless I'm mistaken, the only variant of GFF that explicitly permits Fasta entries is GFF3, and the data you've shown above is definitely not valid GFF3. The mRNA entry seems to correctly use the `=` symbol to separate keys from values in each attribute, but the remaining columns all use spaces and quotation marks (GTF style) instead. This is likely to cause all kinds of problems with software packages that take GFF data as input.

But to more directly address your question, there is no "right" way to connect feature entries to sequence entries in a GFF3 file. There are some *conventions*, but these are all very specific to a small set of tools and community of practice, and it is rare that these conventions are objectively enforced. Some options include:

* Placing genomic nucleotide sequences in the `##FASTA` section, and then using feature coordinates to pull out mature mRNA sequences and translate them *in silico* to derive the corresponding peptide sequence. In this case, the Fasta entries would be genomic sequence identifiers that correspond to the first field/column in each feature entry.
* Placing mature mRNA sequences in the `##FASTA` section, and then directly translating this sequence to derive the corresponding peptide sequence. In this case, you'll want the sequence IDs in the `##FASTA` section to match the `ID` attribute of the corresponding mRNA feature.
* Placing the peptide sequences directly in the `##FASTA` section. In this case, you would either want to have `protein` features declared in the GFF3 file (uncommon) or have some sort of `protein_id` attribute associated with each mRNA feature, in which case you'd want the sequence IDs in the `##FASTA` section to match this `protein_id` attribute.

There are probably other ways of doing it, too. Which way is "right"? The GFF3 specification has nothing to say about it. It depends on the tools you are using to do the subsequent processing and analysis. If you're lucky, those tools might even document what they expect in terms of input.

---

Let me just say, as a post-script, that I think placing sequences and feature annotations in the same file is a horrible idea. I've never liked this about GFF3.

# Answer
> 0 votes
To answer your specific question, you should use the identifier from the first column of the GFF in the ##Fasta section. As you note, this won't work if you are trying to add multiple protein sequence for each seqid. Try using the scaffold DNA sequences in this section and add amino acid translations as an attribute of the annotations.

---
Tags: gff3
---


# Question
Title: Can the canu assembler output a fastq file of the final assembly just like HGAP4?
I have assembled some genome from Sequel PacBio data both with HGAP4 on the SMRT Link interface and using canu on the command line. The HGAP4 assembler outputs a fastq file of the final assembly such that I can see how reliable the sequence is along the genome. On the other hand, canu does no such thing, which means I do not know how much I can trust each nucleotide in the assembly? Is there a way canu can output the fastq file as well or at least some depth information at each base?

# Answer
> 2 votes
I've created an issue on the Canu github repository for this. I'm not aware of any existing functionality to output FASTQ files, but think that this would be a useful feature to have.

It would be possible to create something like this by aligning the original reads to the genome, generating a pileup from that alignment, then determining the sequence likelihood based on that pileup. Unfortunately, I don't have any scripts that I've written to do that... yet.

---
Tags: genome, assembly, fastq, genome-sequencing, canu
---


# Question
Title: Single-cell RNA sequencing (scRNA-seq): filtering cells by transcript counts, how to choose cutoffs?
I am running a notebook with example for the MAGIC algorithm.

In the data preprocessing step, a filtering operation is required to filter out cells with a small count of transcripts.

My question is twofold:

* Why filtering out cells with a small count of transcripts?
* How to decide the cutoff, i.e. the minimum number of transcripts each cell should have to not be removed from the dataset?

I attach some plots from the original notebook, which display some distributions of the single-cell rnaseq dataset.

# Answer
By reading this tutorial:

> Wells with few reads/molecules are likely to have been broken or failed to capture a cell, and should thus be removed.

This answers the first part of the question.

For the second part, i.e. how to decide cutoff, it still seems a bit arbitrary to me, but some authors trim the cells on the left tail of the highest mode, in order to remove other modes (likely due to noise) on the left tail, considering the distribution of cells dependent on the sum of reads (molecules per cell, first plot).
> 1 votes

---
Tags: rna-seq, filtering, single-cell
---


title: Cross-reference with PDB database

text:
# Question
Title: Cross-reference with PDB database
I have a list of several thousand proteins and their UNIPROT IDs. I'm looking for an efficient method of cross-referencing it against the PDB tertiary structure database, and get a list of those proteins with a tertiary structure in the PDB database.

I've tried to BLASTP the list of UNIPROT IDs against the PDB database, using the NCBI BLAST portal but encountered too many errors of "Error: Failed to read the Blast query: Sequence ID not found", making the process of manual filtering not convenient and not efficient.

Any ideas?

Thank you!

# Answer
The UniProt website offers an ID mapping tool that tends to work pretty well:

Then, hit "submit" and you will get the list of corresponding PDB IDs:
> 2 votes

# Answer
You can use the SIFTS mapping: http://www.ebi.ac.uk/pdbe/docs/sifts/

From an article:

> The Structure Integration with Function, Taxonomy and Sequences resource (SIFTS; http://pdbe.org/sifts) is a close collaboration between the Protein Data Bank in Europe (PDBe) and UniProt. The two teams have developed a semi-automated process for maintaining up-to-date cross-reference information to UniProt entries, for all protein chains in the PDB entries present in the UniProt database. This process is carried out for every weekly PDB release and the information is stored in the SIFTS database.

It has no web interface for querying. You just download CSV (or TSV) file with all the mapping data.
> 2 votes

---
Tags: proteins, protein-structure, pdb, uniprot
---


# Question
Title: What are doublets in single cell RNA-seq data?
I am reading The Tabula Muris Consortium et al. (pp).

> In some organs, cells with more than 2 million reads were also excluded as a conservative measure to avoid doublets.

How exactly is a “doublet” defined? For example, is doublet a set of cells sequenced as a single cell (and so, the number of transcripts is double)?

# Answer
> 7 votes
> Is doublet a set of cells sequenced as a single cell?

Yes. Depending on the method of single cell sequencing it may be more or less likely for groups of cells to be captured and barcoded with the same "unique" barcode. This is more likely in split-pool RNA sequencing (e.g. SPLiT-seq), and less likely in cell-capture RNA sequencing (e.g. Fluidigm C1). Doublets can also be created through physical / experimental processes (e.g. from tissues that were not completely dissociated).

Bear in mind that detecting doublets is not as simple as counting for doubling of transcript expression, because the expression profiles are different in different cells. That's why single-cell sequencing is useful in the first place.

# Answer
> 7 votes
"Doublet" is commonly used to describe a droplet in droplet-based sequencing that has captured atleast 2 cells. 10x states their doublet rate to be 0.8% per 1000 cells:

There is a tradeoff between targeted cell capture and doublet rate and lab protocols are usually optimized for a certain expected doublet rate.

One should expect to see close to doubling of transcripts but it is not trivial to infer this computationally.

---
Tags: scrnaseq, single-cell, quality-control
---


title: Lisp in bioinformatics

text:
# Question
Title: Lisp in bioinformatics
After many years of programming in Python, C/C++ and other mainstream languages, I am looking for a new language of pleasure, mainly as a part of mental hygiene. Since I really admire the ideas behind Lisp, including the original paper from the 1950's, I would like to learn some Lisp dialect.

Which Lisp dialects are currently suitable for bioinformatics (e.g., for developing simple command-line tools)? I have the following requirements:

1. Portability across main Unix platforms and reproducibility (this seems to be a huge issue in many Lisp-like languages)
2. Resulting programs being installable using Bioconda (it should be technically possible at least in future if not now)
3. Negligible starting time (this disqualifies Clojure, which requires approximately 10 seconds)

# Answer
I like the idea but I’m not aware of a strong Lisp presence anywhere in bioinformatics. I think a good starting point would be the paper *How the strengths of Lisp-family languages facilitate building complex and flexible bioinformatics applications* by Khomtchouk *&al* (and the Hacker News discussion).

Excluding Clojure from consideration will unfortunately close off a great many modern resources (although I fundamentally agree with your technical reason for doing so). Apart from this, Scheme (in particular R7RS) is probably the most-used Lisp dialect, ahead of Common Lisp. Or at least the most actively developed one. As Pierre mentioned, there’s also GNU Guile but, to be honest, despite it being the default for GNU software tools I don’t think it’s frequently used outside GDB.

That said, Clojure’s startup time is *not* anywhere near 10 seconds on my system. It’s not exactly instantaneous but not much worse than other dynamic languages, despite being backed by the rather sluggish JVM.
> 3 votes

# Answer
> Which LISP dialects are currently suitable for bioinformatics (e.g., for developing simple command-line tools)? I have the following requirements:

I have recently pushed a PR for samtools (mostly for fun) https://github.com/samtools/samtools/pull/752 that uses the **GNU Guile library** https://www.gnu.org/software/guile/

```
./samtools view   -g \
   '(define read-filter (lambda (R) (and (> (hts-read-length R ) 10 ) (integer? (string-contains  (hts-read-seq R) "GCCTTGGCTCTTGTTCC")))))' input.bam

```
> 1 votes

---
Tags: reproducibility
---


# Question
Title: Unable to Install PGDSpider using Ubuntu and wine
I'm running wine v 2.0.3 using ubuntu v. 16.04.2 on a Mac. I'm trying to install a program called PGDSpider2 from here. Essentially the program is supposed to convert different file formats used for bioinformatics.

I'm running it in a 32-bit wine server and I have the Java V.9 and the Microsoft.net packages downloaded on that same server. I've been getting this error.

I have tried the installation on a clean bottle but get the same following error, any thoughts?

```
parallels@ubuntu:~/.wine32/drive_c/Program Files/PGDSpider_2.1.1.3$ WINEPREFIX=~/.wine32 wine PGDSpider2.exe
fixme:heap:RtlSetHeapInformation (nil) 1 (nil) 0 stub
fixme:process:SetProcessDEPPolicy (1): stub
fixme:heap:RtlSetHeapInformation (nil) 1 (nil) 0 stub
fixme:process:SetProcessShutdownParameters (00000380, 00000000): partial stub.
fixme:process:GetProcessWorkingSetSize (0xffffffff,0x73e94c,0x73e950): stub
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e6fc,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
fixme:win:RegisterDeviceNotificationA (hwnd=0x12dd30, filter=0x93e698,flags=0x00000001) returns a fake device notification handle!
fixme:win:UnregisterDeviceNotification (handle=0xcafecafe), STUB!
parallels@ubuntu:~/.wine32/drive_c/Program Files/PGDSpider_2.1.1.3$ fixme:win:EnumDisplayDevicesW ((null),0,0x39edcf0,0x00000000), stub!
fixme:win:EnumDisplayDevicesW ((null),1,0x39edcf0,0x00000000), stub!
X Error of failed request:  BadMatch (invalid parameter attributes)
  Major opcode of failed request:  156 (GLX)
  Minor opcode of failed request:  5 (X_GLXMakeCurrent)
  Serial number of failed request:  241
  Current serial number in output stream:  241
k
fixme:advapi:RegisterEventSourceW ((null),L".NET Runtime Optimization Service"): stub
fixme:advapi:DeregisterEventSource (0xcafe4242) stub

```

# Answer
Don't use Wine at all! You seem to be running an Ubuntu installation inside a MacOS system using parallels. And then a Windows emulation layer (wine) on top of that. This is a very complex system.

However, the software you are trying to install is written in Java and can run perfectly well on any of the operating systems you have available. So there is absolutely no reason to even try to run it via Wine. Just run it natively on either your MacOS or your Ubuntu. As explained in the page you linked to:

Execute PGDSpider GUI:

> * Windows:
> execute the file `PGDSpider2.exe` to start the program
> * Linux:
> execute the command `./PGDSpider2.sh` to start the program
> * Mac and others:
> execute the command `java -Xmx1024m -Xms512m -jar
>  PGDSpider2.jar` to start the program

So ignore the Windows instructions and pick either one of the others.
> 1 votes

---
Tags: format-conversion, linux, phylogenetics, molecular-genetics, mtdna
---