# Question
Title: Mutations in sequences relative to reference (histogram)
I have a reference sequence of amino acids that I compare other sequences to using `Bio.pairwise2.align.globalmx` of Biopython package (one to many compare - one by one \[is there a better way?\])

It's works fine but only returns strings of the alignment with very little metadata.

Now I want to extract the indices of the amino acids that are different than the reference sequence and plot it in some kind of histogram according to protein features:

For example:

* pos 1-36: 3 mutations.
* pos 37-202: 4 mutations.
* pos 203: 1 mutation
* pos 204-660: 0 mutations

Is there any builtin package to do so?

How can I get the data of diff between 2 sequences in python? I mean which AA changed to which AA and at each position. `Bio.pairwise2.align.globalmx` returns some data structures of strings (see below) that can be printed but does not contain the locations in a convenient way.

```
LLGHKVNPFKALYNRIDMK-QVEALVEASKEE
||||||||||||||||||||||||||||||||
LLGHKVNPFKALYNRIDM-RQVEALVEASKEE

```

# Answer
> 2 votes
You can use AlignBuddy to process alignments in GenBank format, which allows you to keep the metadata. From the command line:

```
$: alignbuddy input.gb --generate_alignment mafft | \
   alignbuddy --extract_feature_sequences 'PROT_FEAT_1' | \
   alignbuddy --percent_id

```

If you want to do a large one-to-many comparison and loop through a bunch of features, I'd recommend importing AlignBuddy into a custom python script and calling the relevant functions directly. You can then collect the relevant data and build a histogram with matplotlib.

# Answer
> 1 votes
In answer to the comments, you can do multiple pair wise alignments and then concatenate them to visualise, if MSA is not the way to go. I have a script you can find here which will draw a graph of the Shannon entropy (uncertainty) in a column. If you create a gapped FASTA alignment for each of your pairs, and concatenate it in to a multifasta, the script should still work with it. Though it may be noisy - best just to try it.

JalView might be an option for a one-stop GUI tool, but I’m not sure if it’ll require an MSA instead of multiple pair wise.

---
Tags: sequence-alignment, python, sequence-analysis
---


# Question
Title: Using signal peptide and the expression levels of signal recognition particle in secretome analysis
I have not found any work which investigates assessment of differences in levels of secreted proteins by taking advantage of differential expression of the genes which mediate the secretory pathway.

For example, suppose the gene whose product binds specifically to a certain signal peptide is down-regulated. I would assume that this implies a down-regulation in all the levels of secreted proteins to which it binds.

Is such an analysis possible to date? If so, where?

Thanks in advance.

# Answer
> 1 votes
This is something that Qiagen's Ingenuity Pathway Analysis is meant to solve. Unfortunately, the process of solving pathway dependency problems like these is difficult, and frequently leads to the discovery of super-linked proteins that are seen everywhere (and consequently fairly uninteresting from a biological perspective). The challenge is in filtering out the useless noise of "statistical significance", and identifying the biologically-significant signal. That's not something that computers are particularly good at doing at the moment.

My own suggestion is that you would be better off finding someone who has domain knowledge in the area of secreted proteins, and trying to present results to them in such a way that it makes it easy to interpret the results and ask further questions. If you are that person, perhaps there's a more specific question that you could ask to further investigate the problem that you have.

---
Tags: proteins, differential-expression, networks, interactions
---


# Question
Title: How to measure or assign hydrophobicity score values to individual amino acids of a PDB structure?
I want measure the hydrophobicity of each amino acid within a PDB structure file. Since I have the PDB file I want to consider the 3D information, rather than sequence-only measures such as GRAVY.

I am not interested in measuring the SASA (solvent accessible surface area) or RD (residue depth) - although correlated, it is not the information I require.

How can I comparatively measure the hydrophobicity of the amino acids within a PDB structure file? Is there an algorithm, library or program to do so?

# Answer
Dan is correct but it doesn't answer the question. After much thought, I have finally figured out how to do it, only thanks to Dan's answer.

The problem is that a subsequence is one-dimensional and does not therefore include any of the additional information which is extractable from the three-dimensional structure (when available). Obviously, having the three-dimensional information is much more informative and descriptive for building classifiers, predicting interactions, interaction stability, or other purposes.

1. Find all intramolecular contacts within a radius of N. for the Amino Acid subsequence of interest (such as an Interface). I set N to, e.g. 4 Angstroms.
2. Calculate the descriptive statistics based on the set of hydrophobicity values for all intramolecular contact atoms or amino acids (depending on the granularity). e.g. Make an average of all of the hydrophobicity scores from the selected scale. (I use GRAVY, although there are many others).

To improve this further, the distance value of N can be set more intelligently than as a static value. For example, to the same as the cavity radius.
> 2 votes

# Answer
No, this doesn't makes sense (I think).

The hydrophobicity of an AA is the hydrophobicity of the AA, in whatever context, just like it's molecular weight. It's weight is it's weight, in a structure or in isolation.

I know this seems physically unrealistic, but so is an abstract measure of 'hydrophobicity', it's just the best we can do based on partition coefficients.
> 2 votes

---
Tags: software-recommendation, pdb, 3d-structure, data-mining, computation
---


# Question
Title: LD analysis in PLINK based on reference and a SNP list
I'm new to PLINK and genetics, and getting confused with two PLINK commands for LD analysis:

1. `plink --bfile hapmap --r2 --ld-window-r2 --ld-snp-list --ld-window`
2. `plink --noweb --bfile --clump hapmap --clump inputfile --clump-field P --clump-p1 --clump-p2 --clump-r2 --clump-kb`

It seems to me that for #1, the input files are a reference data and a SNP list file, while for #2, it is a reference data and a .assoc file.

In the output for #2, p-values are included with other results (according to PLINK manual) which is not the case with the output from #1. If I'm interested in finding proxy SNPs for the SNPS on my list in hapmap data, will it be sufficient for me to use the output from #1?

After invoking command #1 using PLINK 1.9, its output gave me 6 columns: "Chromosome number", "chromosomal position", "SNP name" for query/target SNP and 1 column of "R2". I cannot carry out command #1 because I do not have the .assoc file and only have a list of SNPs.

My main question is - what is the difference between these two commands and which of the two should give me the best result for finding proxy SNPs for my list of SNPs.

# Answer
> 2 votes
It looks like the first command is calculating pairwise linkage disequilibrium r2 between your SNPs.

The second command is clumping your dataset based on a p-value over a file of association test results, which I guess is a pre-cursor to computing polygenic risk scores.

Clumping is the process of thinning your SNP list to only keep uncorrelated SNPs so that each representative region has only one marker. In this case you are asking PLINK to keep the marker with the highest p-value.

I'm not sure what your definition of proxy SNPs is but if you are looking for tag SNPs to a SNP of your interest, you should go with command 1 and choose a proper LD r2 threshold.

---
Tags: snp, plink, molecular-genetics
---


# Question
Title: Intersecting two different files with one "master" file based on different columns
I have the following sets of data:

## file1:

```
1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    16855676    16855696    CTGTCCAGCAGAGGGCGGTG

```

## file2:

```
1   15776231    2   5008    G:5002  A:6
1   15776239    3   5008    C:3358  A:14    G:1636

```

## file3:

```
1   16855677    2   5008    A:5003  C:5
1   16855689    3   5008    G:4957  A:41    T:10

```

Most lines have 6 columns but a few lines also have 7 columns which have data that I want to preserve.

Can a script be written which does the following:

if `$2` of file2 comes in between `$2` and `$3` of file1 (between the interval of `$2` and `$3`) and `$2` of file3 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,$7` from file2 and `$5,$6,$7` from file1 and `$2,$3,$4,$5,$6,$7` from file 3, 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

1   15776220    15776240    GTGACCAGCAGGTGTCTCTG    15776239    3   5008    C:3358  A:14    G:1636  16855676    16855696    CTGTCCAGCAGAGGGCGGTG    1   16855689    3   5008    G:4957  A:41    T:10

```

Something like bedtools does but I am not able to use bedtools one this since I don't have two columns in file2 and file3. The intersection is based on $2 in both the files, and all the other columns after it are just copied

(There could be many lines from file2 and file3 intersecting in the same line in file1,so that line could be repeated based on the number of intersections)

# Answer
An awk solution

```
awk '
    FILENAME == ARGV[1] {f2[$2] = $2 OFS $3 OFS $4 OFS $5 OFS $6 OFS $7; next}
    FILENAME == ARGV[2] {f3[$2] = $2 OFS $3 OFS $4 OFS $5 OFS $6 OFS $7; next}
    {
        for (key2 in f2) {
            for (key3 in f3) {
                if (($2 <= key2 && key2 <= $3) && ($5 <= key3 && key3 <= $6)) {
                    print $1,$2,$3,$4, f2[key2], $5,$6,$7, f3[key3]
                }
            }
        }
    }
' file2 file3 file1

```

outputs

```
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776231 2 5008 G:5002 A:6  16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855677 2 5008 A:5003 C:5
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776231 2 5008 G:5002 A:6  16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855689 3 5008 G:4957 A:41 T:10
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776239 3 5008 C:3358 A:14 G:1636 16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855677 2 5008 A:5003 C:5
1 15776220 15776240 GTGACCAGCAGGTGTCTCTG 15776239 3 5008 C:3358 A:14 G:1636 16855676 16855696 CTGTCCAGCAGAGGGCGGTG 16855689 3 5008 G:4957 A:41 T:10

```

For your sample data, all combinations of records in file2 and file3 meet the criteria for the file1 record.
> 4 votes

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


# Question
Title: Inspection of gene expression in scRNA-seq data
I am running the data preprocessing pipeline for scRNA-seq data presented here.

> **3.8.6.1 Gene expression**
>
> In addition to removing cells with poor quality, it is usually a good idea to exclude genes where we suspect that technical artefacts may have skewed the results. Moreover, inspection of the gene expression profiles may provide insights about how the experimental procedures could be improved.
>
> It is often instructive to consider the number of reads consumed by the top 50 expressed genes.

By analyzing my data I have the following most expressed 50 genes (evaluating the total number of transcripts across cells):

**Original**

**Log-transformed**

Since I'll be working on log-transformed data, I assumed that the distributions of the top 50 genes were not excessively skewed. However, as you can see in the first plot, there are some samples (cells) with very high number of transcripts for the top gene. Is it something I should worry about? What other insights can I get from these two figures?

## Edit

Original analysis (figure attached) was:

> The distributions are relatively flat indicating (but not guaranteeing!) good coverage of the full transcriptome of these cells. However, there are several spike-ins in the top 15 genes which suggests a greater dilution of the spike-ins may be preferrable if the experiment is to be repeated.

# Answer
> 3 votes
In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples *prior* to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the point the course instructors are trying to get across with that message box.

Eyeballing the results, you should also exclude any ERCC genes as well, as they are synthetic constructs and not from the target samples (although record their counts, as it may be informative for QC/normalisation purposes); they'll have odd distributions because of this.

ERCC sequences are typically added based on the concentration of DNA in a sample, rather than the number of cells. However with single-cell data (i.e. with known cell counts), it'd be possible to do a proper ERCC-based normalisation, assuming an equal amount were added for each cell. I suspect that knowing a bit more about the ERCC addition protocol would help in understanding why their distribution is so different.

---
Tags: gene, scrnaseq, homework, data-preprocessing
---


# Question
Title: how to count the number of each allele of CNVs?
I have a bunch of files (vcf format) each containing site and sample annotations of a copy number variant (following example is part of one of my files, I am showing 11 individuals but they are nearly 900 individuals)

```
16  32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3

```

I want to count each allele of CNV in each individual(sample) if there are 0, 1 or 2 copies! this is my desired output (rows are individuals):

```
CN0 CN2 CN3 CN4 CN5
0   0   2   0   0
1   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   2   0   0
0   0   1   0   0
0   0   1   0   0
0   0   2   0   0
0   0   0   0   0
0   0   2   0   0
1   0   1   0   0

```

Any idea how to do that?

# Answer
Here's a pure python solution to process your example. It should be straightforward to loop over all the rows of your vcf files

```
my_line='16 32328771    esv3638416;esv3638417;esv3638418;esv3638419;esv3638420  A   <CN0>,<CN2>,<CN3>,<CN4>,<CN5>   100 PASS    AC=27,84,3220,20,26;AF=0.00539137,0.0167732,0.642971,0.00399361,0.00519169;AN=5008;CS=DUP_gs;END=32370968;NS=2504;SVTYPE=CNV;DP=26381;EAS_AF=0.0099,0.0188,0.5427,0.004,0.0089;AMR_AF=0.0101,0.0058,0.4063,0,0.0014;AFR_AF=0.0015,0.0371,0.8374,0.003,0.003;EUR_AF=0.003,0.0089,0.5189,0.005,0.005;SAS_AF=0.0051,0.0031,0.7791,0.0072,0.0072;VT=SV  GT  3|3 3|1 0|3 3|3 3|3 3|0 0|3 3|3 0|0 3|3 1|3 '
my_line=my_line.strip() # remove trailing whitespace
my_fields=my_line.split()
names=my_fields[4].split(',') # extract CNV names
names=[x[1:-1] for x in names] # remove '<>' symbols
print('\t'.join(names))
for sample in my_fields[9:]:
    genotypes=[0]*(len(names)+1) # initialize array of zero counts, including reference
    for allele in sample.split('|'):
        genotypes[int(allele)]+=1
    print('\t'.join([str(x) for x in genotypes[1:]])) # print all columns except reference

```
> 1 votes

---
Tags: python, perl, cnv
---


# Question
Title: Installing RnBeads via bioconductor - .onLoad failed in loadNamespace(), call: NULL
I am trying to install `RnBeads` from bioconductor but the installation of its dependancy `TxDb.Hsapiens.UCSC.hg19.knownGene` fails. The error is:

```
biocLite("RnBeads")
BioC_mirror: https://bioconductor.org
Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.3 (2017-11-30).
Installing package(s) ‘RnBeads’
also installing the dependencies ‘FDb.InfiniumMethylation.hg19’, ‘methylumi’

trying URL 'https://bioconductor.org/packages/3.6/data/annotation/src/contrib/FDb.InfiniumMethylation.hg19_2.2.0.tar.gz'
Content type 'application/x-gzip' length 38996322 bytes (37.2 MB)
==================================================
downloaded 37.2 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/methylumi_2.24.1.tar.gz'
Content type 'application/x-gzip' length 6178909 bytes (5.9 MB)
==================================================
downloaded 5.9 MB

trying URL 'https://bioconductor.org/packages/3.6/bioc/src/contrib/RnBeads_1.10.2.tar.gz'
Content type 'application/x-gzip' length 6437974 bytes (6.1 MB)
==================================================
downloaded 6.1 MB

* installing *source* package ‘FDb.InfiniumMethylation.hg19’ ...
** R
** data
** demo
** inst
** preparing package for lazy loading
Error: package or namespace load failed for ‘TxDb.Hsapiens.UCSC.hg19.knownGene’:
 .onLoad failed in loadNamespace() for 'TxDb.Hsapiens.UCSC.hg19.knownGene', details:
  call: NULL
  error: evaluation nested too deeply: infinite recursion / options(expressions=)?
Error : package ‘TxDb.Hsapiens.UCSC.hg19.knownGene’ could not be loaded
ERROR: lazy loading failed for package ‘FDb.InfiniumMethylation.hg19’

The downloaded source packages are in
    ‘/tmp/RtmpY8Yf4l/downloaded_packages’
installation path not writeable, unable to update packages: doMC, doSNOW,
  foreach, iterators, snow, boot, class, cluster, codetools, foreign,
  KernSmooth, lattice, MASS, Matrix, mgcv, nlme, nnet, rpart, spatial, survival

```

As suggested by the error message, I tried setting `options(expressions=500000)`, which unfortunately did not change the error message. An important note may be that I am under network mount and don't have administrative access rights.

**Edit**

I don't have a .Rprofile file.

Note about `installation path not writeable`: Typing .libPaths() reveals that my home R folder is at the first position (and it is writable). The other paths (`"/usr/local/lib/R/site-library"` , `"/usr/lib/R/site-library"` and `"/usr/lib/R/library"`) are not writable without sudo.

I also tried to delete my home R folder and to reinstall it. The reinstallation of all packages worked for all packages except for the three listed above (same error message). I also tried `biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene", lib.loc=.libPaths()[1], lib=.libPaths()[1])` to specifically address my home R directory, but the same error is showing up (but without the `installation path not writeable` part).

Here is my sessionInfo():

```
R version 3.4.3 (2017-11-30)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Ubuntu 14.04.4 LTS

Matrix products: default
BLAS: /usr/lib/openblas-base/libopenblas.so.0
LAPACK: /usr/lib/lapack/liblapack.so.3.0

locale:
 [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C
 [3] LC_TIME=en_US.UTF-8        LC_COLLATE=en_US.UTF-8
 [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8
 [7] LC_PAPER=en_US.UTF-8       LC_NAME=C
 [9] LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] BiocInstaller_1.28.0

loaded via a namespace (and not attached):
[1] compiler_3.4.3 tools_3.4.3

```

Thank you very much!

# Answer
I have reported the error on support.bioconductor.org but we were not able to solve the problem. It seems that some dependencies are making problems, but as reinstalling my home R folder did not resolve the issue, I have come up with the following workaround:

I downloaded the latest R version from my favorite CRAN mirror, opened an R session with `pathToNewR/bin/R` from my terminal and installed `RnBeads` there. This works fine.

Note: I had to make sure that `pathToNewR/library` is the only entry in `.libPaths()`, as per default, my home R folder was also contained in `.libPaths()`.

As a consequence, I am independent from the (maybe broken) R version and packages on my server and can finally work with `RnBeads` :)
> 2 votes

---
Tags: r, bioconductor
---


# Question
Title: Difference between de novo transcriptome assembly methods
I have been looking around (including read the original papers) to understand what is essentially the difference between StringTie in non-reference based mode (de novo) and Trinity de novo assembly. I understand that in the genome-guided nature of StringTie, we are supported by the genome information (i.e. taking into account the mapped reads). In addition to that, we also have an option to either run using the gene annotation file or not.

Simply put, am I correct if I say that Trinity is purely looking at the overlapping reads, then assemble them, whereas StringTie (non-reference based) is looking at the mapped reads that are either close or overlapping to each other? Is there any other thing I'm missing here (I understand it's not that simple, but I'm trying to make it as simple/intuitive as possible)?

Also if that's the case, would it be safe to say that generally non-reference based StringTie is preferred over Trinity and reference based StringTie (as gene annotation can lead to some "bias", hence preventing discovery of novel TSS, exon length or even make unnecessary duplicates with small differences)?

# Answer
Your intuition is correct. stringTie is just looking at clumps of alignments and how they might relate to each other (either due to spliced alignments or proximity). Trinity is doing the more computationally difficult task of finding parts of reads that overlap with each other and trying to link them together into longer sequences.

Whether it makes sense to use something like stringTie or not depends mostly on the quality of the reference sequence. If you're working on an organism with a high-quality reference sequence (pretty much anything with a lot of papers using it) then there's no reason to bother with Trinity. If your reference genome additionally has a high-quality annotation (human, mouse, etc.) then you should use it, as anything "novel" that the program is likely to output otherwise is probably just an artefact. If you either don't have an annotation or it's known to be highly problematic then it makes sense to omit it when using stringTie.

Without seeing evidence to the contrary, I don't buy the argument that using a good annotation is actually biasing the output in a problematic way. For well-studied organisms, anything reported as "novel" these days likely doesn't exist, or is otherwise so fleeting that it's biologically irrelevant. The arguments of bias sound like something that was come up with once someone didn't get as many "novel" results as they were hoping.
> 3 votes

# Answer
It's really a misnomer to call StringTie's non-reference based mode 'de-novo.' It's still using the reference genome sequence to guide the transcript assembly, it's just not using the reference annotation. Trinity is truly de-novo in that it it assembles the transcript from the overlap of the reads without mapping them to a reference genome sequence.

If you look at the precision numbers from the StringTie paper, you'll see that it maxes out around 80% under ideal circumstances with simulated data. This means that if you're working with a well annotated genome where most transcripts have already been discovered, the vast majority of the 'novel' transcripts you find will likely turn out to be artifacts.
> 3 votes

---
Tags: rna-seq, sequence-annotation, assembly, transcriptome
---


# Question
Title: Why doesn't Biopython AlignIO.read() recognise the 'mauve' format?
On this Biopython tutorial, they describe how to import a multiple sequence alignment in the Mauve (XMFA: extensible multi fasta format). So I imported the AlignIO module:

```
from Bio import AlignIO
alignment = AlignIO.read(open("alignment.xmfa"), "mauve")

```

but I get the following error:

```
ValueError
Traceback (most recent call last) <ipython-input-30-396e54791d7a> in <module>()
----> 1 alignment = AlignIO.read(open("/Users/cr517/Documents/alignment.xmfa"),
"mauve")

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in read(handle, format, seq_count, alphabet)
    422     iterator = parse(handle, format, seq_count, alphabet)
    423     try:
--> 424         first = next(iterator)
    425     except StopIteration:
    426         first = None

/Users/cr517/anaconda/lib/python3.6/site-packages/Bio/AlignIO/__init__.py
in parse(handle, format, seq_count, alphabet)
    366                                                 seq_count=seq_count)
    367         else:
--> 368             raise ValueError("Unknown format '%s'" % format)
    369
    370         # This imposes some overhead... wait until we drop Python 2.4 to fix it

ValueError: Unknown format 'mauve'

```

# Answer
You're using version 1.68 or older. Mauve support was added in 1.70.
> 5 votes

---
Tags: sequence-alignment, python, file-formats, errors
---


# Question
Title: How to filter ribosomal RNA from scRNA-seq data
I want to filter out ribosomal RNA from scRNA-seq data (downloaded from here). Is there a list of known ribosomal RNA?

The only solution I found is SortMeRNA, however it works with raw sequencing data afaik, while I already have a matrix with transcript counts for each cell. I searched for a comprehensive list of rRNAs but I didn't find any.

# Answer
> 3 votes
The rRNA genes in that dataset are `Rn45s` and `Rn4.5s`.

BTW, you have gene counts, not transcript counts.

# Answer
> 6 votes
In the paper mentioned, we used the `ScaleData` function in Seurat to regress out the number of reads, Rn45s abundance, and percent ribosomal gene transcripts. Ribosomal genes were found with the regular expression `^Rp[sl][[:digit:]]`.

```
tiss <- ScaleData(object = tiss, vars.to.regress = c("nReads", "percent.ribo","Rn45s"))

```

Here's a fuller notebook, and we'll have a better organized repository soon.

---
Tags: rna-seq, scrnaseq, rna, ribosomal
---


# Question
Title: Hba-a2 has two copies on UCSC genome browser
I was browsering the UCSC viewer of the genome at this region.

UCSC genes label two genes as Hba-a2 while GENCODE thinks the left one is Hba-a1

Which annotation is correct?

# Answer
Gencode is correct (when in doubt, always assume Gencode/Ensembl is correct when it differs from UCSC). What presumably happened is that the UCSC annotation pipeline noticed that Hba-a1 and Hba-a2 are nearly identical (they produce the same protein and have minimal nucleotide difference) and just misannotated.
> 2 votes

---
Tags: ucsc, genome-browser
---


title: How to show only my read on IGV?

text:
# Question
Title: How to show only my read on IGV?
I have a large BAM file on my IGV. I'd like to visualize a read alignment. I have the read name, and I don't care anything else in the file.

IGV has an option "Select by name" (see screenshot). The problem is that IGV doesn't hide the alignments not matching the name. I'd still have to **manually navigate** on IGV with my mouse to find my read.

That was too cumbersome. What I'm doing now is make a new alignment BAM file. For example,

`samtools view -h L.V.51.bam | grep '@\|HWI-D00119:228:CB08RANXX:8:2116:4095:45403' | samtools view -bS - | samtools sort - > new.bam`

I'd load the new alignment file on my IGV.

If I have 20 reads, I'd have to create 20 new BAM files... Is there anything I can do on IGV that I'm not aware of?

# Answer
Partially a shameless plug... I've written a genome browser, ASCIIGenome, entirely controlled through the command line. As such, it makes relatively easy to deal with questions like yours. For example, you could do:

```
ASCIIGenome bam/WW00292a.bam  # Load bam file(s)

```

Filter for reads matching the regex passed to `-i` option:

```
[h] for help: grep -i '72789|57688|70785|56299|27116|12086'

```

(Optionally, show sam records as strings)

```
[h] for help: print

```

Example result:

For command help type `grep -h` and for more filtering options see also `awk -h`.

Hope this helps and feel free to post questions and open issues on github!
> 7 votes

---
Tags: igv
---


# Question
Title: Error in .checkedCall : subset indices out of range
After normalization the way it is described here:

https://bioconductor.org/packages/3.7/bioc/vignettes/scran/inst/doc/scran.html

Using `quickCluster()` function and running `computeSumFactors()` I am trying to run `calculateQCMetrics()` the way it is described here:

https://f1000research.com/articles/5-2122/v2

```
sce <- calculateQCMetrics(sce, feature_controls=list(Mt=is.mito))

```

And I am getting an error:

> Error in .checkedCall(cxx\_colsum\_subset, exprs\_mat, is\_feature\_control) : subset indices out of range

Does anybody know how to correct it? What is going wrong here?

# Answer
> 2 votes
Just posted the question and found out why the error was happening. I had two datasets: raw, and filtered one with half of the rows as raw. I did normalization on the filtered one. Then I calculated `is.mito` using the unfiltered, raw dataset. So, it produced indexes that were outside of the range of the dataset.

---
Tags: r, scrnaseq, scran
---


# Question
Title: Drawbacks of upper quartile normalization for scRNA-seq data
I would like to use Upper Quartile normalization for scRNA-seq data defined as:

> The upperquartile (UQ) was proposed by (Bullard et al. 2010). Here each column is divided by the 75% quantile of the counts for each library. Often the calculated quantile is scaled by the median across cells to keep the absolute level of expression relatively consistent. A drawback to this method is that for low-depth scRNASeq experiments the large number of undetected genes may result in the 75% quantile being zero (or close to it). This limitation can be overcome by generalizing the idea and using a higher quantile (eg. the 99% quantile is the default in scater) or by excluding zeros prior to calculating the 75% quantile.

In particular, I'd use the 99% quantile to avoid normalizing by small number or zero. Is there any other drawback for using it on scRNA-seq data?

* Bullard, James H, Elizabeth Purdom, Kasper D Hansen, and Sandrine Dudoit. 2010. “Evaluation of Statistical Methods for Normalization and Differential Expression in mRNA-Seq Experiments.” BMC Bioinformatics 11 (1). Springer Nature: 94. doi:10.1186/1471-2105-11-94.

# Answer
My primary concern with using the top ~1% or so in upper quantile normalization is that it's going to be prone to the same robustness issues that RPKM/FPKMs run in to. That is, if for whatever technical reason you have to have a fair bit of variability in a couple very highly expressed genes (typically rRNAs, but one can imagine other genes) and the set of genes going into the UQ normalization is composed primarily of those then the normalization results will simply reflect irrelevant technical variability.

Of course, this all assumes that it's even appropriate to normalize scRNAseq data with a standard method like this. In some cases it is (e.g., when comparing control vs. treatment with a largely homogeneous cell population), but if you're dealing with very different cell types then I worry that the standard method is going to wash away things like differences in bulk RNA content between cell types. Presuming there are spike-ins, they should have the same normalization applied and then checked to ensure there's nothing greatly amiss. Having said that, if you're going to be doing tSNE and want to look for cell types, then odds are good you'll still be able to see most of them (after all, it's unlikely that RNA quantity is the only difference between them, though not removing that difference would make things easier).

# The problem of global shifts in expression

The concept of global shifts in expression is not unique to scRNAseq, in fact the most common example of this is Myc up/down regulation in relation to cancer. To explain the issue, let's suppose that you have the following (very small) matrix of counts:

```
       sample1 sample2
gene1        1       2
gene2       10      20
gene3       15      32
gene4      100     198

```

As you can see, sample2 has roughly twice as many reads for each gene as sample1, so any normalization is going to cut all of its counts in half (or the reverse for sample1, or perhaps a mixture of the two...depending on the method). That's usually fine, since the assumption with all common normalization methods is that differences like this are simply due to differences in sequencing depth between the samples. But what happens if sample2 actually has twice as much RNA as sample1? Well, then you've just normalized away all of the differences between the samples (I encourage readers to go through the literature on Myc, since that's a very nice example of this, where a mostly-global shift in expression led to apparent down-regulation of genes with constant expression). Such a situation doesn't occur very commonly in bulk RNAseq, but when it does you need to use something like spike-ins (dosed to match cell number or DNA quantity) to maintain the global difference.

If this is only rarely an issue with bulk RNAseq, then why is it more likely to be an issue in scRNAseq? The answer is that bulk RNAseq benefits from the averaging out of many different cell-types, such that ensemble averages/medians don't tend to change much due to treatments...but this isn't going to be the case at the single-cell level, since cells can vary vastly in size (compare a sperm and an oocyte...every RNA should be differentially expressed between them). This difference in base RNA content between different cell types will be washed away by common normalization methods, which can lead to either losing differences or creating differences in genes where there shouldn't be (again, see examples of this in the Myc literature). Whether this really becomes an issue will depend on what your goals are with your data and how much of a global expression difference you actually have and the type of analysis you want to perform.

At the end of the day, just remember that common RNAseq normalization methods assume that (A) there's no global expression change between samples and that (B) if there are very many DE genes that the direction of change is not heavily skewed in one direction. If either of these assumptions are violated then you can run into issues.
> 4 votes

---
Tags: normalization, scrnaseq
---


# Question
Title: getting KeyError on certain residues in BioPython: ATOM vs. HETATM lines
When I want to load atoms for residuum `(' ', 85, ' ')` in chain *A* for protein **3tmm** with `structure[0][chain_id][residue]`, but I'm getting `KeyError` exception:

```
Traceback (most recent call last):
  [...]
  File "/home/pok/Checkouts/bilkoviny/loader.py", line 118, in _get_residue_atoms
    return list(self.structure[0][chain_id][residue])
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Chain.py", line 74, in __getitem__
    return Entity.__getitem__(self, id)
  File "/usr/lib64/python2.7/site-packages/Bio/PDB/Entity.py", line 40, in __getitem__
    return self.child_dict[id]
KeyError: (' ', 85, ' ')

```

I assume this is because this residuum do not have "ATOM" lines in PDB file:

```
[...]
ATOM    704  CD1 ILE A  84      19.348  -4.839  18.138  1.00 47.08           C
ANISOU  704  CD1 ILE A  84     5862   6108   5917   -216  -1513    609       C
HETATM  705  N   MSE A  85      21.619  -0.130  16.236  1.00 31.89           N
ANISOU  705  N   MSE A  85     5143   3072   3900   -766  -1011    962       N
HETATM  706  CA  MSE A  85      22.277   0.675  15.226  1.00 40.76           C
ANISOU  706  CA  MSE A  85     6708   3869   4911   -972   -729   1162       C
HETATM  707  C   MSE A  85      23.146   1.786  15.818  1.00 44.98           C
ANISOU  707  C   MSE A  85     7191   4141   5757  -1268   -558   1109       C
HETATM  708  O   MSE A  85      24.257   2.045  15.338  1.00 43.69           O
ANISOU  708  O   MSE A  85     7055   3855   5692  -1611   -151   1153       O
HETATM  709  CB  MSE A  85      21.244   1.265  14.269  1.00 46.50           C
ANISOU  709  CB  MSE A  85     8043   4309   5316   -709   -948   1345       C
HETATM  710  CG  MSE A  85      21.849   1.954  13.052  1.00 54.03           C
ANISOU  710  CG  MSE A  85     9604   4929   5995   -907   -623   1638       C
HETATM  711 SE   MSE A  85      20.503   2.755  11.877  1.00 84.35          SE
ANISOU  711 SE   MSE A  85    14348   8336   9366   -486  -1038   1884      SE
HETATM  712  CE  MSE A  85      21.690   3.538  10.518  1.00 83.79           C
ANISOU  712  CE  MSE A  85    15063   7865   8907   -915   -414   2326       C
ATOM    713  N   SER A  86      22.636   2.453  16.846  1.00 43.76           N
ANISOU  713  N   SER A  86     5234   6599   4794  -1778   -946    195       N
[...]

```

I have been told "MSE is same as MET, it is just with selenium instead of sulphur" and that I probably need to load these atoms (I'm approximating secondary structure element position, so I want individual atom coordinates).

Is there some way how to get atoms for this residue in this protein via BioPython?

# Answer
You can find the answer in The Biopython Structural Bioinformatics FAQ:

> **What is a residue id?**
>
> This is a bit more complicated, due to the clumsy PDB format. A residue id is a tuple with three elements:
>
> * The **hetero-flag**: this is 'H\_' plus the name of the hetero-residue (e.g. 'H\_GLC' in the case of a glucose molecule), or 'W' in the case of a water molecule.
> * The **sequence identifier** in the chain, e.g. 100
> * The **insertion code**, e.g. 'A'. The insertion code is sometimes used to preserve a certain desirable residue numbering scheme. A Ser 80 insertion mutant (inserted e.g. between a Thr 80 and an Asn 81 residue) could e.g. have sequence identifiers and insertion codes as follows: Thr 80 A, Ser 80 B, Asn 81. In this way the residue numbering scheme stays in tune with that of the wild type structure.
>
> The id of the above glucose residue would thus be ('H\_GLC', 100, 'A').

So in your case it would be `('H_MSE', 85, ' ')`.
> 1 votes

---
Tags: python, pdb
---


# Question
Title: Approximate matches with bedtools intersect
I am using `bedtools` to test whether a set of variant calls produced by a new method is missing any variants from a set of known validated variants.

```
bedtools intersect -v -a newcalls.vcf -b goldstandard.vcf

```

It's working as intended for SNV calls, but it occasionally misses indels due to a single base pair discrepancy. It's usually an alignment artifact that doesn't affect the interpretation of the variant.

My first thought was to make a BED file with intervals that extend a few bp beyond each variant call, and comparing *this* BED file to the known variants. However, I was wondering whether approximate matching is possible with `bedtools intersect`?

# Answer
> 2 votes
You can use `bedtools intersect` but first use `bcftools norm`. bcftools norm will left-align and normalize indels to get both vcf files consistent with each other. This should resolve the matching if they are truly matching variants.

# Answer
> -1 votes
You could use `bedtools closest` with `-d` option (*-d: In addition to the closest feature in B, report its distance to A as an extra column.*). Indels that have no match in file A would be those where the closest indel in file B is more than *x* bases away.

---
Tags: vcf, bedtools
---


title: GSEA- parsing error for .cls file

text:
# Question
Title: GSEA- parsing error for .cls file
I am creating a .cls file for my GSEA analysis and see the below error. Can anyone suggest a solution?

I create the file using TeXworks. The file looks like this:

```
37 2 1
#CD69neg CD69pos
CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

```

## Error Details

```
---- Full Error Message ----
There were errors: ERROR(S) #:1 Parsing trouble edu.mit.broad.genome.parsers.Par ...

---- Stack Trace ----

of exceptions: 1
------Number of items found in cls data 10 is not equal to the number specified on the header line 37------

edu.mit.broad.genome.parsers.ParserException:
Number of items found in cls data 10 is not equal to the number specified on the header line 37 at
edu.mit.broad.genome.parsers.ClsParser.doSanityChecks(ClsParser.java:388) at
edu.mit.broad.genome.parsers.ClsParser._parse_genecluster_style_categorical(ClsParser.java:252) at
edu.mit.broad.genome.parsers.ClsParser.parse(ClsParser.java:226) at
edu.mit.broad.genome.parsers.ParserFactory._readTemplates(ParserFactory.java:342) at
edu.mit.broad.genome.parsers.ParserFactory.readTemplate(ParserFactory.java:293) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:753) at
edu.mit.broad.genome.parsers.ParserFactory.read(ParserFactory.java:726) at
edu.mit.broad.genome.parsers.ParserWorker.doInBackground(ParserWorker.java:52) at
javax.swing.SwingWorker$1.call(SwingWorker.java:295) at
java.util.concurrent.FutureTask.run(FutureTask.java:266) at
javax.swing.SwingWorker.run(SwingWorker.java:334) at
java.util.concurrent.ThreadPoolExecutor.runWorker(ThreadPoolExecutor.java:1149) at
java.util.concurrent.ThreadPoolExecutor$Worker.run(ThreadPoolExecutor.java:624) at java.lang.Thread.run(Thread.java:748)

```

# Answer
You seem to have forgotten the space after the #. It should be as in here (Just for future readers know where to find the specs). Also it should be a single line (Maybe this error was introduced while pasting here or it has been a duplicate because there are 39 samples in the 3 line and they should be 37).

Either adding 2 more samples to the initial count and removing the newline of the third line:

```
39 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

```

Or removing the third line:

```
37 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

```
> 3 votes

---
Tags: gsea
---


# Question
Title: Is there anything similar to GSEA for locus-based (instead of of gene-based) data?
As the question states, I am interested in an analysis similar to Gene Set Enrichment Analysis (ranked gene sets) but focused on locus-level data instead of genes.

To explain in greater detail: I have a set of genomic coordinates from DNA methylation data that are scored based on their contribution to components from a non-negative matrix factorization. I am interested in better understanding the biological underpinnings of what these NMF components represent. I know I can use Locus Overlap Analysis tools like GREAT or LOLA to tell me what GO/Functions are enriched in a particular BED file, but this type of analysis doesn't take into account the scores/rank of the loci in the BED file as would be done in an analysis akin to GSEA.

Does anyone know of any tools/methods for analyzing functional enrichment while taking advantage of feature/locus ranks/scores with locus-level data?

I know I probably didn't explain this very clearly, so please let me know if I can clarify anything for you. Thanks

# Answer
I tinkered with a program GSEA-SNP quite a few years back, which claims that it does a similar ranking procedure with SNPs. It carries out its procedure by first linking SNPs with genes, then running an algorithm similar to GSEA. There's a bit more detail in the paper.

Unfortunately this is a space I've got a bit further away from in recent years, so I don't know if anyone has carried out the legwork to generalise the algorithm for sets of any type where each item within each set fits a continuous (but not necessarily normal) distribution.
> 2 votes

# Answer
You can use GSEA-like methods for any type of data. The basis of the test is to analyze if there is a preference for a group (a gene set) in the ranking of the other variable.

So you first need to find where the functions you are interested are in the genome (To associate a range/position with a given function). Then using this association you can test if the ranking of the locus is associated with a specific function that is in a given area.

I haven't seen anything like that implemented but you could do it yourself using biomaRt to retrieve data from where a gene (or a function) is located and fgsea to test which functions are associated with your score.
> 0 votes

---
Tags: gsea, go-enrichment, functional-annotation
---


# Question
Title: Parse out exon coordinates from bed file for each gene
This is my bed file for all the exon coordinates I want to take out all the exons of a given gene. Let say I have gene in chr 1 which starts from `chr1 11868 12227` so I want to parse out all the exons that comes in between `11868 12227`.

This is my small subset:

```
cat exon.bed | head -10
chr1    11868   12227   +   exon
chr1    11871   12227   +   exon
chr1    11873   12227   +   exon
chr1    12009   12057   +   exon
chr1    12178   12227   +   exon
chr1    12594   12721   +   exon
chr1    12612   12697   +   exon
chr1    12612   12721   +   exon
chr1    12612   12721   +   exon
chr1    12974   13052   +   exon

```

How do I parse out?

I use mostly R and bit of shell script but I'm not sure if I can use R, may a few lines of perl or shell script can help me solve my problem.

# Answer
There are multiple ways to go about it. On the command line you can make a 1 line BED file:

```
chr1   11868   12227

```

And then `bedtools intersect` with it.

In R, you could load your original BED file and use GenomicRanges:

```
library("GenomicRanges")
bed = read.delim("foo.bed", header=F)  # Rename this
# N.B., BED files use 0-based coordinates, I've switched to 1-based here
exons = GRanges(seqnames=levels(bed$V1),
                ranges=IRanges(start=bed$V2+1, end=bed$V3),
                strand=bed$V4)
# Again, I switched to 1-based coordinates...you don't have to do that
region = GRanges(seqnames="chr1",
                 ranges=IRanges(start=11869, end=12227),
                 strand="+")
theExonsYouWant = subsetByOverlaps(exons, region, type="within")

```

You'll need to either specify a strand in your region or use `ignore.strand=T` in `subsetByOverlaps()`. I hope you don't have any genes that overlap on the same strand, since they'll prove problematic.
> 5 votes

# Answer
Via BEDOPS `bedops -n` and Unix I/O streams:

```
$ echo -e "chr1\t11868\t12227" | bedops -n 1 exon.bed - > answer.bed

```

Or, if you have your genes in a BED file called `genes.bed`:

```
$ bedops -n 1 exon.bed genes.bed > answer.bed

```

If you have your genes in some other format, like GFF or GTF, you can use `gff2bed` or `gtf2bed`, e.g.:

```
$ bedops -n 1 exon.bed <(gff2bed < genes.gff) > answer.bed

```

Or:

```
$ bedops -n 1 exon.bed <(gtf2bed < genes.gtf) > answer.bed

```

The file `answer.bed` will contain exons that do not overlap a gene annotation.
> 3 votes

---
Tags: r, bed, bedtools, perl
---


# Question
Title: Is this a low complexity region in our human genome?
I have a screenshot where many of my reads were aligned to a region that I **suspect** a low complexity region. Although you can't see, all those reads are **clipped** in the cigar strings. Sample cigar strings: *64S32M29S, 74S32M18S* etc... Consequently, the actual sequence mapped to the genome is **less** than the read length.

I have a feeling that my alignments are bogus because of the complexity and clipping, but I'm not sure if the complexity is indeed low. All I see is a few "T" between a bunch of "A". How to define a low complexity region?

**Q:** Is this a low complexity region in the human hg38 genome? Would that be fair if I report "*the alignments are likely be errors due to the low complexity and alignment clipping*"?

# Answer
Yes, that's a short low complexity region wedged between a SINE and an snRNA. More importantly, your alignments have MAPQ of 0 (that's why they're filled with white in IGV), which will happen if they map equally well to multiple locations. Without looking at the sequence one can use that alone to determine that these are not trustworthy mappings.
> 7 votes

# Answer
Yes, this is a low complexity region.

Regions are considered low complexity (or having a simple sequence) when they contain an abundance of a single base, or an abundance of short tandem repeats. The simplest is a tandem repeat of a single base (e.g. `AAAAAAAAAAAAAAAA`), also called a homopolymer. In your case there is a 3-unit tandem repeat of `TAAAAA`, and surrounding As (with a few intermingled Ts). Depending on the genetic history of this region, it might be that those other As and Ts are derived from the same `TAAAAA` unit, with some slipping that has occurred during replication to create a cryptic (i.e. non-perfect) repeat structure.
> 6 votes

---
Tags: genome
---


# Question
Title: Improve scRNA-seq dataset for further analysis
I got a dataset from `C.Elegans scRNA-seq` paper:

`GSM2599701_Gene.count.matrix.celegans.cell.Rdata` in `GSE98561_RAW.tar`

The dataset is `40 000 x 68 000`, where rows represent genes and columns - cells. So, I took it and tried to process myself to build an `scRNA-seq` pipeline. Here is what I did:

1. I filtered out the genes that have zero counts in all of the cells and the dataset was reduced to `29 000 x 68 000`
2. I removed all of the cells with counts `< 100` in all of the genes - the dataset became `29 000 x 66 000`
3. Then, because the dataset was too big to run normalization even on the cluster with `120Gb` RAM (because there are multiple distinct types of cells, first clustering needs to be done), I selected just even columns and ran normalization with `29 000 x 33 000` dataset (`UMI_count`):

```
library(scran)
library(scater)
sce <- newSCESet(countData=UMI_count)
clusters <- quickCluster(sce)
sce <- computeSumFactors(sce, clusters=clusters, positive=TRUE)

```

After running the code above I decided to check whether the data is fine and so I ran:

```
> summary(sizeFactors(sce))
Min.    1st Qu.  Median  Mean   3rd Qu.    Max.
0.0000  0.0000  0.0000  0.0717  0.0000 33.3900

```

I also ran `PCA` on the normalized dataset and it looks like that:

It seems to me that the normalized dataset is terrible and I need to do some further processing before I could run further analysis. What else could I do to improve it? How to filter it? There are no `spike-ins` and maybe `200` mitochondrial genes. The approach described here does not work, probably because the majority of the cells have low number of genes expressed:

I tried removing low-abundance genes after normalization, but it seems that most of them are going to be removed:

```
>ave.counts <- rowMeans(counts(sce))
>keep <- ave.counts >= 1
>sum(keep)
109

```

Should I filter all of the cells - columns - for up to `500` total gene expression count instead of `100`? Is it a good idea? I can't think of anything else.

# Answer
> 3 votes
I usually use a minimum of 200 genes/cell for relaxed filters (500 if possible). For genes I use presence in at least 10% of the cells.

A few other filters you can try are:

* High percentage of reads from only top50 genes. This means you haven't captured a lot of genes and the diversity of the cell won't be useful.
* High percentage of mitochondrial genes. This has been reported to be linked to bad quality cells. The idea is that mRNA is depleted because of early cell death (or other outer membrane injury) leading to having only mito mRNA left.
* Minimum UMI per cells. Similar to minimum genes but more dependent of the protocol used.
* UMI per gene/cell. This can be interesting to detect a potential bias from amplification and is also dependent on the protocol.

The best would be to compare all your samples using violin plots with those different filters and choose your thresholds from there. A quick example (using Seurat):

Pre filtering:

Post filtering:

# Answer
> 5 votes
I don't think you can conclude that the dataset is terrible based on that PCA. Depending on the specific protocol, each scRNA-seq dataset is going to be very different. Unlike bulk RNA-seq where all the samples are going to be of very similar quality, individual cells will be highly variable. For example, the most basic QC metric is the number of reads. If some of your bulk RNA-seq samples have 10X more reads than other samples, that's a cause for concern. In scRNA-seq experiments, that's not at all surprising. The PCA exposes the major sources of variation. In scRNA-seq experiment, those are probably going to be technical rather than biological, so the output may not make much sense. Furthermore, even if all the cells are of perfect quality, PCA is not really optimal for sparse data with many dropouts like you see with scRNA-seq.

Regardless of PCA, those library sizes seem very small. Based on your plot, there are barely any cells with over 10K reads. I understand this is C elegans, but that seems somewhat low. That may be why you are having trouble with normalization. It would help with understanding the data if you changed the X-axis on those plots so that you can actually see where the peak is. Right now, it looks like 0, but that can't be the case.

# Answer
> 1 votes
Given that this is single-cell sequencing results, you should be using `rowSums` instead of `rowMeans` for filtering genes. Most genes won't be expressed at a single cell level, but there should be a lot of genes that have some expression. The zero counts are probably better represented as "missing" rather than zero.

Unfortunately that feature makes running PCAs difficult, because PCAs don't work with missing values. The missing values may need to be imputed in order to generate a better PCA.

disclaimer: I don't \[typically\] work with single cell data, so don't have a good handle on the methods for looking at it.

---
Tags: r, scrnaseq, scran
---


title: Find number of possibilities

text:
# Question
Title: Find number of possibilities
We have a column which has following format : (1, 1, 1, 0, 1).

Here is the definition of what means two columns are **compatible**.

> Using the notation Oi to denote the collection of rows possessing a 1 in their i-th element, we conclude that one of three possibilities must be true: Oi ⊆ Oj , Oj ⊆ Oi , or Oi and Oj are disjoint (the first two cases include the possibility that Oi = Oj). If columns i and j satisfy this condition, then we call them compatible

The question is : **How many possible columns of length 5 are compatible with the column (1, 1, 1, 0, 1)?**

You can read more here where problem is described in page 14-16.

Thanks in advance !

# Answer
> 2 votes
Let $n$ be the number of rows (in your case $n=5$) and $|O\_j|$ the cardinality of $O\_i$.

We can calculate the number of compatible columns as sum of the number of $O\_i \subseteq O\_j$, $O\_j \subseteq O\_i$ and $O\_i, O\_j$ disjoint.

number of $O\_i \subseteq O\_j$: Here we need to find the number of supersets of our $O\_i$. That means we need to find the number of sets which contain the indices in $O\_i$ and for each index not in $O\_i$, we have two choices: It is contained, or it is not contained in $O\_j$. Therefore the number of supersets is $2^{n-|O\_i|}$, which is 2 to the power of elements not in our set.

number of $O\_i, O\_j$ disjoint: Here all elements of $O\_i$ are fixed not to be in $O\_j$ and for each index not in $O\_i$, it can again be in $O\_j$ or not. That means that the number of possible sets is again $2^{n-|O\_i|}$.

number of $O\_j \subseteq O\_i$: Here we need to find the number of subsets of $O\_i$ and this is the number of elements in its power set, which is $2^{|O\_i|}$.

So finally, this means that the number of compatible columns can be calculated as: $2^{n-|O\_i|} + 2^{n-|O\_i|} + 2^{|O\_i|} - 2$.

Why -2? We have counted 2 sets twice:

* The empty set, which is a subset of $O\_i$ and also counted in the number of disjoint sets
* and $O\_i$, which we counted in $O\_i \subseteq O\_j$ and in $O\_j \subseteq O\_i$

This means for your examle, we have $16 + 2 + 2 - 2 = 18$ compatible sets.

**Note:** The all zero vector (which corresponds to $O\_i=\\{\\}$) is a special case for which the formula above does not work, as all disjunct sets and all supersets are counted twice. For this special case, the number of compatible columns is simply $2^n$.

---
Tags: gene, genome, homework
---


# Question
Title: Creating a new reference genome B by changing genome A with mismatches from BAM file without long reads
This question is slightly related to this one: Improve a reference genome with sequencing data

Only in my case, I have a starting genome reference A, of a short phage genome, and reads from a related genome B, for which I haven't been able to find a reference.

I took the Illumina 2x75bp paired end reads from genome B and ran Velvet on it (I ran Ray, but it didn't complete). It gives me a fragmented collection of contigs. The contigs align 99-98% to reference genome A. I don't have long reads.

My question is: with the short reads and the contigs that I have, can I transform the reference genome A into a reference genome B?

PS: no, I don't have long reads, and there is no budget for them. The budget for long reads is zero and will remain zero. Even if we could sequence long reads for free, we would still not do it because of the lab work.

Tools that I have tried so far:

ococo:

```
$HOME/ococo/ococo -i contigs.bwamem.bam -c 1 -f $HOME/Genomes/T2_phage/starting.fasta -F consensus.fasta --vcf-cons - -V variants.vcf

```

It produced a consensus.fasta, and checking how much it changes from the reference, it changes. There are one or more changes in about 33% of the 100bp windows:

```
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep '\ \|\ ' | wc -l
1272
diff -y <(seqtk seq -l 100 $HOME/Genomes/T2_phage/starting.fasta) <(seqtk seq -l 100 consensus.fasta) | grep -v '\ \|\ ' | wc -l
419

```

# Answer
Even though ococo does not seem to handle big indels well, it's a good enough tool for a first try.
> 0 votes

# Answer
You should be able to use something like Pilon to convert from one genome to another, assuming the changes are all fairly small local changes (but that's almost never the case).

I find your resistance to long-read data interesting. It's important to consider the cost of bioinformatics as well. Sequencing a phage genome with long reads will take less than an hour, and probably produce a few reads that cover the entire genome. With existing short-read data, those full-length reads should be able to be fixed up within a day.

I spent about half a year trying to fix a eukaryotic genome assembly that was assembled only from short reads (using Pilon, and L-RNA-scaffolder), and I was only able to reduce the number of contigs down by about half (and even at that, I wasn't convinced that contigs were being properly joined). With long reads and another couple of weeks of computer work, the contigs dropped down to about a tenth of their original number.

If you can use a pipette and have about half an hour spare, you can do rapid long-read sequencing on the MinION.
> 2 votes

---
Tags: sam, reference-genome, snp
---


title: Heteroplasmy Calculation

text:
# Question
Title: Heteroplasmy Calculation
For calculating heteroplasmy in mtDNA I was using MitoSeek, the tool provides position wise heteroplasmy percentage. So, will it be wise to add up all sites to get a position-wise percentage in a single sample, or calculate it by adding up all counts (variants) from an aligned file and divide it by total read count?

# Answer
Heteroplasmy at *any* site (as long as it's a reliable variant) suggests that there is mitochondrial heteroplasmy, it doesn't matter where it is. In that case, the "mitochondrial heteroplasmy" statistic would simply be the maximum value at any \[reliable\] location.

Due to the low amount of variation in the mitochondrial genome, and because it's a haploid genome, I don't think that a calculation like "average heteroplasmy" makes much sense for mitochondria.
> 1 votes

---
Tags: mtdna, heteroplasmy
---


# Question
Title: Filter out outliers of the scRNA-seq (heterogenous cells)
I am new to data science. I have a dataset of single-cell gene expression from multiple cell types in *C. Elegans*. The dataset is from the paper Comprehensive single-cell transcriptional profiling of a multicellular organism

My main question is, **Which approaches should I use for filtering out bad cells in this case when we have multiple cell types in the dataset?**

So far I tried to filter out genes that have too high mitochondrial genes content following the Bioconductor “simpleSingleCell” workflow.

However, the tutorial specifically says that the method of filtering out based on mitochondrial genes most probably will not work when the dataset has multiple cell types:

> Analyzing all cell types together would unnecessarily inflate the MAD and compromise the removal of low-quality cells, at best; or lead to the entire loss of one cell type, at worst.

Any suggestions would be greatly appreciated.

# Answer
> 2 votes
From what I known, there is no clear consensus in the field and it depends on the type of cells you are interrogating.

However, if mitoRNA/endogenousRNA ratio does not fit for your purposes, other option is to **check at the total number of genes/transcripts detected** in each cell. In this way, you can filter out cells that have considerable less genes/transcripts detected than the rest since this may be indicative of bad quality cells for any reason (apoptosis, RNA degradation, sequencing itself, etc). For instance, to **put a threshold** of \>1000 genes/transcripts in a cell to be considered for further analysis.

Consider also to take a look at the number of mapped reads for each cell, as cells with low mapped reads respect to the rest could be potentially problematic.

---
Tags: rna-seq, filtering
---


title: What is the ICGC normalized_read_count?

text:
# Question
Title: What is the ICGC normalized_read_count?
I downloaded gene expression data (`exp_seq`) from the ICGC file browser.

For each sample and gene, the file contains a `normalized_read_count`.

What is that value? I couldn't find any information on the ICGC website. The values are definitly too low for TPM.

# Answer
> 5 votes
By reading this thread on seqanswers and by comparing the data to TCGA, I figured out

* `raw_read_count` is the read count which you use as input for e.g. `DESeq2`. It has been estimated using RSEM
* `normalized_read_count` is equivalent to the `scaled_estimate` from TCGA. This is the *estimated fraction of transcripts made up by a given gene*, as estimated by RSEM. Multiplying this value with 1e6 yields the TPM.

---
Tags: rna-seq, public-databases
---


# Question
Title: What is the detection limit of protocol in scRNA-seq
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.

I used to think that the detection limit of a protocol was the average maximum number of molecules detected across cells in a scRNA-seq experiment. However, by working on the dataset provided by The Tabula Muris Consortium et al. (pp) I noticed that RNA raw counts often spikes to more than 10 molecules.

What is the actual meaning of detection limit of a protocol in scRNA-seq?

# Answer
> 1 votes
The detection limit refers to the ability to detect a given gene/transcript. In this case, they're saying that a given transcript requires 5-10 copies in order to be likely detected (in other words, the loss rate is 80-90%, which is directly related to a previous question of yours).

---
Tags: scrnaseq
---


# Question
Title: Does the MAPQ=0 fraction of a BAM file depend on the insert sizes?
When doing Illumina 2x150bp sequencing of genomic DNA, and after aligning the reads to GRCh38, does the percentage of the non-N fraction of the human genome as MAPQ=0 depend on the insert sizes of the genomic fragments?

This is, for two identical samples with identical final coverage, sample A having an average insert size of 250bp and sample B having an average of 450bp, would the fraction of MAPQ=0 change between the two?

# Answer
Yes, as a general rule of thumb mappability increases with insert size (up to a limit) and read length. Whether this will actually occur in a given case will depend more on how randomly the sequencing samples from the genome to begin with (i.e., if the library prep happens to select for/against high mappability regions then the insert size won't much matter). The basic reason behind this is that if one end of a read aligns to a repetitive element that the likelihood of being able to use the other end as an anchor will increase a bit if it's further away. Now having said that, Illumina sequencers in particular like a certain fragment length for efficient sequencing, so you can't go too crazy on increasing the insert sizes and still get decent output (unless you modify your library prep).
> 4 votes

# Answer
> This is, for two identical samples with identical final coverage, sample A having an average insert size of 250bp and sample B having an average of 450bp, would the fraction of MAPQ=0 change between the two?

For human, I *expect* sample B will have fewer MAPQ=0 reads (I don't have concrete data to prove). A large part of human genome consists of ALU repeats. The average size of an ALU is ~320bp. If your insert size is ~450bp, you can map reads into most of isolated ALUs in principle. With ~250bp insert, fragments in the middle of ALUs won't be mappable uniquely.

Note that to see the difference, you need to use a mapper that sufficiently takes the advantage of paired-end information.
> 0 votes

---
Tags: sam, illumina, mapq
---


title: Plot percentage of genome covered

text:
# Question
Title: Plot percentage of genome covered
Given a aligned bam file (wgs/bwa-mem) what steps do i need to perform to generate a plot as seen below (from this paper):

We are trying to see how much of genome was covered using pacbio and illumina for cow.

---

You mean remove this left plot?

# Answer
> 5 votes
The general steps are:

1. Install deepTools
2. Run plotCoverage on your BAM files
3. Remove the left-most plot in photoshop/gimp/imagemagick

# Answer
> 1 votes
Devon gives a practical answer using the great deepTools package. These should be the theoretical, lower level steps behind the plots:

* Compute read depth at each genomic position, using for example `samtools mpileup` with appropriate flags to accept/reject reads. You get a table like:

```
Pos       Depth
chr1:1    0
chr1:2    10
chr1:3    11
...
chrX:100  10

```
* For each unique value *x* in the *Depth* column count how many genomic positions have depth \>= *x*. You can use Unix tools like `sort ... | uniq -c | awk ...` for this. You get a table like:

```
Depth   N bases with depth >= x   % bases with depth >= x
0       (all genomic positions)   100 %
1       1000                      90 %
2       995                       85 %
.
.
.
N       2                         0.1 %

```
* Now just plot column `Depth` (x) vs `% bases with depth >= x` (y) with R or you favourite plotting tool.

Did I get it right?

---
Tags: wgs
---


# Question
Title: What is the meaning of the Hmmer profile?
I want to understand what the hmmer Profile used in hmmer.

consider the profile

```
HMMER2.0  [converted from 3.1b2]
NAME  aRNH_profile
LENG  121
ALPH  Amino
RF    no
CS    no
MAP   yes
NSEQ  34
DATE  Tue Apr 26 21:24:57 2016
XT       -8455     -4  -1000  -1000  -8455     -4  -8455     -4
NULT        -4  -8455
NULE       656  -1722     98    419   -333    475  -1125    239    250    947  -1073   -271    -50   -338    114    451    113    430  -2131   -717
HMM        A      C      D      E      F      G      H      I      K      L      M      N      P      Q      R      S      T      V      W      Y
         m->m   m->i   m->d   i->m   i->i   d->m   d->d   b->m   m->e
           -6      *  -7932
     1    -214   -916   -565    694    862  -1448    512   -838    115  -1042    170   1348  -1481    369   -158   -290   -159   -625    -52   1675     1
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378    -18      *
     2     285    990  -2213  -1678   2305  -1819   -698    219  -1289   -140    545  -1167  -1835   -973  -1322   -424   1025    762   -559    377     2
     -    -210   -336    220    -37   -342    377     73   -669    210   -511   -732    404    423     70     27    439    205   -414   -161   -176
     -     -18  -6890  -7932   -894  -1115   -701  -1378      *      *
[....]

```

what exactly does each component mean? By definition a HMM $\lambda$ consist of an alphabet $V$ in my case twenty amino acids. Different states $S$ in our case main, insertion and deletion (beginning and end). And two probability matrices of state transitions and emission in certain state.

now if I look at the profile I can see that it has the allowed transitions
`m->m m->i m->d i->m i->i d->m d->d b->m m->e` and all the amino acids, but what does each row in the table mean. I checked the manual from the documentation and it says rather ironically

> If you were to look at this file (and you don’t have to – it’s intended for HMMER’s consumption, not yours), you’d see \[part of my output\]

are these number translations of the probabilities in $\[0,1\]$ into another interval $\[-x,+y\]$? what is the meaning of the three rows under `A 1`

```
-124
-210
 -18

```

Is this the transition matrix? It would be great if someone could explain this or point to the correct reference.

# Answer
> 1 votes
The tutorial in the documentation does indeed state what you quoted but a full description of the file format is given in Section 8 of the documentation, starting on page 106.

To summarise, as you noted, the number columns A-K is the alphabet size (the number of amino acids). The number of rows underneath is the number of nodes, which corresponds to the number of match states. The three numbers (-214, -210 and -18 in your example) are: the match emission, insert emission, and state transition values.

Caveat: The HMMER profile HMM format changes frequently between HMMER versions (even between minor versions) so there is no guarantee this information is correct for all the older file formats.

---

Quoting the relevant part of the documentation on page 109:

> **main model section**
>
> All the remaining fields are mandatory. The first two lines in the main model section are atypical.‡ They contain information for the core model’s BEGIN node. This is stored as model node 0, and match state 0 is treated as the BEGIN state. The begin state is mute, so there are no match emission probabilities. The first line is the insert 0 emissions. The second line contains the transitions from the begin state and insert state 0. These seven numbers are: B → M1, B → I0, B → D1; I0 → M1, I0 → I0; then a 0.0 and a ’*’, because by convention, nonexistent transitions from the nonexistent delete state 0 are set to log 1 = 0 and log 0 = −∞ = ‘*’.
>
> The remainder of the model has three lines per node, for M nodes (where M is the number of match states, as given by the LENG line). These three lines are (K is the alphabet size in residues):
>
> **Match emission line**
>
> The first field is the node number (1 . . . M). The parser verifies this number as a consistency check (it expects the nodes to come in order). The next K numbers for match emissions, one per symbol, in alphabetic order. The next field is the MAP annotation for this node. If MAP was yes in the header, then this is an integer, representing the alignment column index for this match state (1..alen); otherwise, this field is ‘-’. The next field is the CONS consensus residue for this node. If CONS was yes in the header, then this is a single character, representing the consensus residue annotation for this match state; otherwise, this field is ‘-’. The next field is the RF annotation for this node. If RF was yes in the header, then this is a single character, representing the reference annotation for this match state; otherwise, this field is ‘-’. The next field is the MM mask value for this node. If MM was yes in the header, then this is a single ’m’ character, indicating that the position was identified as a masked position during model construction; otherwise, this field is ‘-’. The next field is the CS annotation for this node. If CS was yes, then this is a single character, representing the consensus structure at this match state; otherwise this field is ‘-’.
>
> **Insert emission line**
>
> The K fields on this line are the insert emission scores, one per symbol, in alphabetic order.
>
> **State transition line**
>
> The seven fields on this line are the transitions for node k, in the order shown by the transition header line: Mk → Mk+1, Ik, Dk+1; Ik → Mk+1, Ik; Dk → Mk+1, Dk+1. For transitions from the final node M, match state M + 1 is interpreted as the END state E, and there is no delete state M + 1; therefore the final Mk → Dk+1 and Dk → Dk+1 transitions are always * (zero probability), and the final Dk → Mk+1 transition is always 0.0 (probability 1.0).

---
Tags: hidden-markov-models
---


title: PLINK clump behavior on missing SNPs?

text:
# Question
Title: PLINK clump behavior on missing SNPs?
I have a long list of autoimmune-associated SNPs, and I want to boil it down so that I get one SNP representing each LD block. I chose to use PLINK's --clump option for this. I'm roughly following this tutorial (but analyzing my own data).

I don't have the full data from the original GWAS's, so I am using 1000 Genomes data as the reference for LD estimation. But, the 1000 Genomes data don't have all the SNPs I need. I end up a few hundred warnings:

`Warning: 'rs12345678' is missing from the main dataset, and is a top variant.`

How does PLINK handle this situation? More importantly, how should I handle this situation? I don't want to omit those SNPs, but I also don't want to take them all, for the same reason that I am doing this in the first place.

# Answer
Most likely, the GWASs that generated your summary statistics used other imputation panels than 1000G, like HRC. Clearly, PLINK can't estimate the LD for SNPs that aren't found in the reference, and most likely just ignores them. I think the easiest solution is to use whatever individual-level dataset you have to estimate the correlations - even if it's small or from a different population. The estimate will be biased, but you don't really need high precision for clumping.

Then again, you can skip PLINK and try searching for SNP prioritization methods based on summary statistics - it's an area of active development, and something might be available to suit your needs.
> 6 votes

# Answer
You can use imputation for guessing at what the missing SNPs might be based on known LD patterns in populations. This procedure will give you an idea of whether the recombinational history of genomic regions is important, rather than specific SNPs. Given that you're looking at LD-based measures, this should be okay (although the logic behind the calculation is a bit circular).

If you don't have a big compute cluster available, the Michigan Imputation Server is a free imputation service that can be used for academic purposes. Data is fed into the service as bgzip-compressed VCF files (output from `plink` as an uncompressed VCF file, then use `bgzip` to compress the files), and returned as a larger dataset with imputed SNPs derived from the HRC or 1kg reference panels. There's a bit of initial legwork required: you need to make sure that the variants as defined in the VCF file have the same location and strand as defined in the reference panel, otherwise you'll wait weeks in the queue and end up empty handed. To avoid hassle, I'd recommend removing any complementary SNPs (e.g. A/T, C/G), then flipping the strand for any remaining genotypes where the alleles don't match up.

I've recently done this again. Here's a summary of the things required to get MIS imputation working properly from a PLINK dataset:

1. Retrieve site data for HRC panel
2. Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
3. extract out variant names (remove duplicated rs numbers)
4. change 'RS' to 'rs' in plink BIM file
5. filter PLINK file against variants
6. adjust position of variants to fit with HRC panel
7. sort HRC variants by rs number
8. identify variants that need flipping
9. flip variants using PLINK, create combined VCF file
10. create VCF index for bcftools
11. extract per-chromosome VCF files
12. convert chromosome '23' back to chromosome 'X'
13. upload to MIS, remembering to separate out autosomal and X-chromosomal files

And a full example:

```
## Retrieve site data for HRC panel
wget 'ftp://ngs.sanger.ac.uk/production/hrc/HRC.r1-1/HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz'
## Generate list of complementary variants, rename chrX to chr23 to fit PLINK expectation
zgrep -v '^#' HRC.r1-1.GRCh37.wgs.mac5.sites.vcf.gz | cut -f 1-5 | \
  awk '{if(!($4 ~ /[AT]/ && $5 ~ /[AT]/) && !($4 ~ /[CG]/ && $5 ~ /[CG]/)){print}}' | perl -pe 's/^X/23/' > HRC_complementary_variants.tsv
## extract out variant names (remove duplicated rs numbers)
cut -f 3 HRC_complementary_variants.tsv | grep '^rs' | uniq -u > HRC_complementary_variants_namesOnly.txt
## change 'RS' to 'rs' in plink BIM file
perl -i -pe 's/RS/rs/' plinkFileHg19.bim
## filter PLINK file against variants
plink --bfile plinkFileHg19 --extract HRC_complementary_variants_namesOnly.txt --make-bed --out noComp_plinkFileHg19
## adjust position of variants to fit with HRC panel
plink --bfile noComp_plinkFileHg19 --update-chr HRC_complementary_variants.tsv 1 3 \
  --update-map HRC_complementary_variants.tsv 2 3 --make-bed --out repositioned_nc_plinkFileHg19
## sort HRC variants by rs number
sort -k 3,3 HRC_complementary_variants.tsv > rsSorted_HRC_complementary_variants.tsv
## identify variants that need flipping
join -1 3 -2 2 <(pv rsSorted_HRC_complementary_variants.tsv) <(sort -k 2,2 repositioned_nc_plinkFileHg19.bim) | \
  awk '{if(!($4 == $9) && !($4 == $8)){print $1}}' > toFlip_HRC_complementary_variants.txt
## flip variants using PLINK, create combined VCF file
plink --bfile repositioned_nc_plinkFileHg19 --flip toFlip_HRC_complementary_variants.txt --recode vcf bgz --out flipped_rnc_plinkFileHg19
## create VCF index for bcftools
bcftools index flipped_rnc_plinkFileHg19.vcf.gz
## extract per-chromosome VCF files
for x in $(seq 1 23); do echo ${x};
  bcftools filter -r ${x} flipped_rnc_plinkFileHg19.vcf.gz -O z -o chr${x}_flipped_rnc_plinkFileHg19.vcf.gz;
done
## convert chromosome '23' back to chromosome 'X'
zcat chr23_flipped_rnc_plinkFileHg19.vcf.gz | perl -pe 's/ID=23/ID=X/;s/^23\t/X\t/' > chrX_flipped_rnc_plinkFileHg19.vcf
bgzip chrX_flipped_rnc_plinkFileHg19.vcf

```
> 3 votes

---
Tags: snp, gwas, plink
---


# Question
Title: What algorithms exist for reverse protein folding?
Are there any algorithms already in existence for reverse protein folding/docking (e.g. taking a desired surface map/docking profile and generating DNA which will produce a protein to fold into it?)

Protein folding seems straightforward (though computationally intensive) and protein docking is likewise straightforward, but thus far I haven't come across any protein unfolding algorithms, though I might just be searching for the wrong keywords.

Any advice in this area?

# Answer
The researchers behind FoldIt have the occasional *design challenges* that use human game players to help design proteins to fit a structure.

However, even with a known protein sequence there's still a lot of labwork involved in trying to find a good working DNA sequence. Just generating the DNA sequence is a difficult task due to redundancy in the genetic code (i.e. 64 codons for 20 amino acids + stop signals). When that is combined with cell-specific tRNA frequencies that can influence the timing of amino acid incorporation (and therefore how a protein folds as it is created), creating a sequence that reliably fits a shape ends up being quite difficult.
> 1 votes

---
Tags: docking
---


# Question
Title: How can I extract gene names for a metabolic pathway from KEGG?
*Note: this question has also been asked on Biostars*

I need the get the list of gene names involved in glycolysis (to put an example). Not manually, I need to do this in a script. Ideally with Python. I think KEGG is the proper database to do this, but any other (free) suggestion would do.

How can I do it?

# Answer
> 6 votes
I whipped this little script up using the KEGG API:

```
#!/usr/bin/env python3
import urllib.request
import re
import sys
pathway = 'hsa00010' # glycolysis
url = "http://rest.kegg.jp/get/" + pathway
with urllib.request.urlopen(url) as f:
    lines = f.read().decode('utf-8').splitlines()
    want = 0
    for line in lines:
        fields = line.split()
        ## The list of genes starts here
        if fields[0] == 'GENE':
            want = 1
            ## The line with GENE is different
            print(fields[2].rstrip(';'))
        ## We reached the next section of the file
        elif want == 1 and re.match('^\S', line):
            sys.exit();
        ## We're still in the list of genes
        if want == 1 and len(fields)>1:
            print(fields[1].rstrip(';'))

```

# Answer
> 6 votes
using the REST API of togows: one can fetch a json structure of the genes associated to a kegg pathway. eg: http://togows.org/entry/kegg-pathway/hsa00010/genes.json

```
[
  {
    "3101": "HK3; hexokinase 3 [KO:K00844] [EC:2.7.1.1]",
    "3098": "HK1; hexokinase 1 [KO:K00844] [EC:2.7.1.1]",
    "3099": "HK2; hexokinase 2 [KO:K00844] [EC:2.7.1.1]",
    "80201": "HKDC1; hexokinase domain containing 1 [KO:K00844] [EC:2.7.1.1]",
    "2645": "GCK; glucokinase [KO:K12407] [EC:2.7.1.2]",
    "2821": "GPI; glucose-6-phosphate isomerase [KO:K01810] [EC:5.3.1.9]",
    "5213": "PFKM; phosphofructokinase, muscle [KO:K00850] [EC:2.7.1.11]",
    "5214": "PFKP; phosphofructokinase, platelet [KO:K00850] [EC:2.7.1.11]",
    "5211": "PFKL; phosphofructokinase, liver type [KO:K00850] [EC:2.7.1.11]",
    "2203": "FBP1; fructose-bisphosphatase 1 [KO:K03841] [EC:3.1.3.11]",
    "8789": "FBP2; fructose-bisphosphatase 2 [KO:K03841] [EC:3.1.3.11]",
    "230": "ALDOC; aldolase, fructose-bisphosphate C [KO:K01623] [EC:4.1.2.13]",
    "226": "ALDOA; aldolase, fructose-bisphosphate A [KO:K01623] [EC:4.1.2.13]",
    "229": "ALDOB; aldolase, fructose-bisphosphate B [KO:K01623] [EC:4.1.2.13]",
    "7167": "TPI1; triosephosphate isomerase 1 [KO:K01803] [EC:5.3.1.1]",
    "2597": "GAPDH; glyceraldehyde-3-phosphate dehydrogenase [KO:K00134] [EC:1.2.1.12]",
    "26330": "GAPDHS; glyceraldehyde-3-phosphate dehydrogenase, spermatogenic [KO:K10705] [EC:1.2.1.12]",
    "5232": "PGK2; phosphoglycerate kinase 2 [KO:K00927] [EC:2.7.2.3]",
    "5230": "PGK1; phosphoglycerate kinase 1 [KO:K00927] [EC:2.7.2.3]",
    "5223": "PGAM1; phosphoglycerate mutase 1 [KO:K01834] [EC:5.4.2.11]",
    "5224": "PGAM2; phosphoglycerate mutase 2 [KO:K01834] [EC:5.4.2.11]",
    "441531": "PGAM4; phosphoglycerate mutase family member 4 [KO:K01834] [EC:5.4.2.11]",
    "2027": "ENO3; enolase 3 [KO:K01689] [EC:4.2.1.11]",
    "2026": "ENO2; enolase 2 [KO:K01689] [EC:4.2.1.11]",
    "2023": "ENO1; enolase 1 [KO:K01689] [EC:4.2.1.11]",
    "5315": "PKM; pyruvate kinase, muscle [KO:K00873] [EC:2.7.1.40]",
    "5313": "PKLR; pyruvate kinase, liver and RBC [KO:K12406] [EC:2.7.1.40]",
    "5161": "PDHA2; pyruvate dehydrogenase alpha 2 [KO:K00161] [EC:1.2.4.1]",
    "5160": "PDHA1; pyruvate dehydrogenase alpha 1 [KO:K00161] [EC:1.2.4.1]",
    "5162": "PDHB; pyruvate dehydrogenase (lipoamide) beta [KO:K00162] [EC:1.2.4.1]",
    "1737": "DLAT; dihydrolipoamide S-acetyltransferase [KO:K00627] [EC:2.3.1.12]",
    "1738": "DLD; dihydrolipoamide dehydrogenase [KO:K00382] [EC:1.8.1.4]",
    "160287": "LDHAL6A; lactate dehydrogenase A like 6A [KO:K00016] [EC:1.1.1.27]",
    "92483": "LDHAL6B; lactate dehydrogenase A like 6B [KO:K00016] [EC:1.1.1.27]",
    "3939": "LDHA; lactate dehydrogenase A [KO:K00016] [EC:1.1.1.27]",
    "3945": "LDHB; lactate dehydrogenase B [KO:K00016] [EC:1.1.1.27]",
    "3948": "LDHC; lactate dehydrogenase C [KO:K00016] [EC:1.1.1.27]",
    "124": "ADH1A; alcohol dehydrogenase 1A (class I), alpha polypeptide [KO:K13951] [EC:1.1.1.1]",
    "125": "ADH1B; alcohol dehydrogenase 1B (class I), beta polypeptide [KO:K13951] [EC:1.1.1.1]",
    "126": "ADH1C; alcohol dehydrogenase 1C (class I), gamma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "131": "ADH7; alcohol dehydrogenase 7 (class IV), mu or sigma polypeptide [KO:K13951] [EC:1.1.1.1]",
    "127": "ADH4; alcohol dehydrogenase 4 (class II), pi polypeptide [KO:K13980] [EC:1.1.1.1]",
    "128": "ADH5; alcohol dehydrogenase 5 (class III), chi polypeptide [KO:K00121] [EC:1.1.1.1 1.1.1.284]",
    "130": "ADH6; alcohol dehydrogenase 6 (class V) [KO:K13952] [EC:1.1.1.1]",
    "10327": "AKR1A1; aldo-keto reductase family 1 member A1 [KO:K00002] [EC:1.1.1.2]",
    "217": "ALDH2; aldehyde dehydrogenase 2 family (mitochondrial) [KO:K00128] [EC:1.2.1.3]",
    "224": "ALDH3A2; aldehyde dehydrogenase 3 family member A2 [KO:K00128] [EC:1.2.1.3]",
    "219": "ALDH1B1; aldehyde dehydrogenase 1 family member B1 [KO:K00128] [EC:1.2.1.3]",
    "501": "ALDH7A1; aldehyde dehydrogenase 7 family member A1 [KO:K14085] [EC:1.2.1.3 1.2.1.8 1.2.1.31]",
    "223": "ALDH9A1; aldehyde dehydrogenase 9 family member A1 [KO:K00149] [EC:1.2.1.3 1.2.1.47]",
    "221": "ALDH3B1; aldehyde dehydrogenase 3 family member B1 [KO:K00129] [EC:1.2.1.5]",
    "222": "ALDH3B2; aldehyde dehydrogenase 3 family member B2 [KO:K00129] [EC:1.2.1.5]",
    "220": "ALDH1A3; aldehyde dehydrogenase 1 family member A3 [KO:K00129] [EC:1.2.1.5]",
    "218": "ALDH3A1; aldehyde dehydrogenase 3 family member A1 [KO:K00129] [EC:1.2.1.5]",
    "84532": "ACSS1; acyl-CoA synthetase short chain family member 1 [KO:K01895] [EC:6.2.1.1]",
    "55902": "ACSS2; acyl-CoA synthetase short chain family member 2 [KO:K01895] [EC:6.2.1.1]",
    "130589": "GALM; galactose mutarotase [KO:K01785] [EC:5.1.3.3]",
    "5236": "PGM1; phosphoglucomutase 1 [KO:K01835] [EC:5.4.2.2]",
    "55276": "PGM2; phosphoglucomutase 2 [KO:K15779] [EC:5.4.2.7 5.4.2.2]",
    "2538": "G6PC; glucose-6-phosphatase catalytic subunit [KO:K01084] [EC:3.1.3.9]",
    "57818": "G6PC2; glucose-6-phosphatase catalytic subunit 2 [KO:K01084] [EC:3.1.3.9]",
    "92579": "G6PC3; glucose-6-phosphatase catalytic subunit 3 [KO:K01084] [EC:3.1.3.9]",
    "83440": "ADPGK; ADP dependent glucokinase [KO:K08074] [EC:2.7.1.147]",
    "669": "BPGM; bisphosphoglycerate mutase [KO:K01837] [EC:5.4.2.11 5.4.2.4]",
    "9562": "MINPP1; multiple inositol-polyphosphate phosphatase 1 [KO:K03103] [EC:3.1.3.80 3.1.3.62]",
    "5105": "PCK1; phosphoenolpyruvate carboxykinase 1 [KO:K01596] [EC:4.1.1.32]",
    "5106": "PCK2; phosphoenolpyruvate carboxykinase 2, mitochondrial [KO:K01596] [EC:4.1.1.32]"
  }
]

```

# Answer
> 3 votes
Answers from Biostars:

GeneSCF 'prepare\_database' module will extract all the pathways with corresponding genes as simple table format in plain text file. Check 'Two step process' (sub-heading) from the provided link (First step will do your job).

---

You can use KEGG REST API for extracting gene names from the pathway. Python request library can be used to make the REST calls to KEGG REST API and query information. Please follow this tutorial if you are a newbie to the Python programming and REST web services.

---
Tags: python, pathway
---


# Question
Title: Where can I find a list of SNP rs IDs that belong to the X chromosome?
Where can I find a list of SNP rs IDs that belong to the X chromosome?

I'm trying to exclude the SNPs on X-chromosome in a data set, but I don't know where these SNPs come from.

I tried to use dbSNP to search, it returns 26933249 SNPs, but the webpage seems to refuse me to download the information.

# Answer
You just searched for `X[Chromosome]`, you didn't specify a species. Presumably, your dataset comes from a specific species, so you should limit your search to that species. For instance, for human, you can search dbSNP for `X[Chromosome] human[organism]`. Then, click on the `send to` link at the top left and download the data:

This is veeeeeery slow, however. I suggest that you instead download the data for chrX from dbSNP's FTP site (this assumes your species of interest is human): ftp://ftp.ncbi.nih.gov/snp/organisms/human\_9606/chr\_rpts/chr\_X.txt.gz

Once you have this file, you can use it to filter your dataset:

1. Collect all rsIDs from the file:

```
zcat chr_X.txt.gz | awk '/^[0-9]/{print "rs"$1}' > chrX.rsIDs

```
2. Filter your dataset.

You haven't shown the format of your dataset, so if it is just a list of rsIDs like this:

```
$ cat dataset
rs77669712
rs199706910
rs2799940
rs12205639
rs10958205
rs200778574

```

You can remove any that are on chrX with:

```
$ grep -vxFf chrX.rsIDs dataset
rs77669712
rs199706910
rs12205639
rs10958205
rs200778574

```

That will print all rsIDs from the text file `dataset` that are not also in the text file `chrX.rsIDs`. Note the `-x` option for `grep` which tells it to only consider matches that cover the entire line. This ensures that `rs1234` is not removed when you're looking for `rs123`.
> 6 votes

---
Tags: snp, phylogenetics
---


title: bed file to gtf conversion using awk

text:
# Question
Title: bed file to gtf conversion using awk
I have a bed file like this:

```
chr1    12227   12612   LOC102725121    0   +
chr1    12721   13220   LOC102725121    0   +
chr1    12227   12612   DDX11L1 0   +
chr1    12721   13220   DDX11L1 0   +
chr1    14829   14969   WASH7P  0   -
chr1    15038   15795   WASH7P  0   -
chr1    15947   16606   WASH7P  0   -
chr1    16765   16857   WASH7P  0   -
chr1    17055   17232   WASH7P  0   -
chr1    17368   17605   WASH7P  0   -

```

I am trying to convert into gtf file for downstream analysis. Gtf file should look like this:

```
chr1   myIntrons   feature 3073253   3074322   .   +   .     gene_id "intron1"; transcript_id "intron1"
chr1   myIntrons   feature 3075858   3079636   .   +   .     gene_id "intron2"; transcript_id "intron2"

```

For the gene\_id part in the 9th column in the double quotes I want my 4th column and in transcript\_id I want combination of column1,column2,and column3 so that my output looks like this: Showing for the first line

```
chr1    myIntrons   intron  12228   12612   .   +   .   gene_id "LOC102725121"; transcript_id "chr1:12227-12612"

```

I am using this code but under double quotes it is not expanding `$1 $2` etc

```
awk -v OFS="\t" '{print $1,"myIntrons","intron", $2+1, $3, ".", $6, ".", "gene_id $4; transcript_id $1:$2-$3"}' file.bed| head

```

# Answer
> 2 votes
Instead of the `"gene_id ..."` bit, you want `sprintf("gene_id \"%s\"; transcript_id \"%s:%s-%s\"", $4, $1, $2, $3)`.

---
Tags: awk, gtf
---


# Question
Title: ethnicity check either from bam or vcf files
What tool could I use to check the ethnicity of a human bam or vcf file? I would like to use the results as a QC check to know whether a given sample or set of samples match the ethnicity information annotated in the meta data or not.

Seen so far:

The examples work, but starting from a GRCh38 .vcf file (from GATK UG) gives me an error:

```
Read    mytest.vcf   :1 Individual(s)Exception in thread "main" java.lang.NumberFormatException: For input string: ""
        at java.lang.NumberFormatException.forInputString(NumberFormatException.java:65)
        at java.lang.Integer.parseInt(Integer.java:592)
        at java.lang.Integer.valueOf(Integer.java:766)
        at Individual.get_Individual(Individual.java:73)
        at similarity.main(similarity.java:80)

```

It won't work with GATK UG .vcf files straightaway, it needs to reprocess .bam files using the ANGSD pipeline:

This tool does not give ethnicity checks but only genotype checks. Running all-vs-all vcfs vs bams can create a semi-matrix of similarities.

# Answer
The main difficulty here is the use of GRCh38. Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference. It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the best.

Let's suppose you go with overlapping rsIDs, and a list of those rsIDs, one per line is in 'rsids.txt'. (Due to strand flips between reference builds, you may want to restrict this list so that all allele codes match, and it's also reasonable to throw out A/T and C/G SNPs.) Then, the following process would work:

1. Download ADMIXTURE (https://www.genetics.ucla.edu/software/admixture/download.html ), plink (https://www.cog-genomics.org/plink/1.9/ ), and/or chr1-chr22 of 1000 Genomes phase 3 (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/20130502/ ) if you don't already have them.
2. Use plink to extract the overlapping higher-MAF rsIDs from each 1000 Genomes phase 3 VCF ("plink --vcf ... --extract rsids.txt --maf 0.05 --make-bed --out ..."), and then merge the resulting per-chromosome filesets ("plink --merge-list ... --out merged\_phase3\_subset").
3. ADMIXTURE prefers a dataset with around 100k variants. Unless you have a small list of overlapping rsIDs (in which case you had better use another method), you should have more than that. plink's LD-pruning function is a good way to select a subset for ADMIXTURE's use: "plink --bfile merged\_phase3\_subset --indep-pairwise 500kb 1 0.2; plink --bfile merged\_phase3\_subset --extract plink.prune.in --make-bed --out admixture\_data"

Adjust the 0.2 threshold as necessary to keep the right number of variants.

4. Run ADMIXTURE in unsupervised mode ("admixture admixture\_data.bed 5 -j8"; adjust the -j parameter depending on the number of processor cores). This generates an admixture\_data.5.P file with population allele frequencies, and an admixture\_data.5.Q file with sample population assignments. Verify that the .Q file corresponds to the 1000 Genomes phase 3 'superpopulations'.
5. Convert your data to plink-format if necessary, keeping just the overlapping rsIDs and 'downgrading' to GRCh37 coordinates. Make sure your variants are sorted in GRCh37-coordinate order, and the allele order also matches that in your reference dataset ("--a2-allele admixture\_data.bim 6 2" during your final --make-bed operation will do the trick); otherwise ADMIXTURE won't do the right thing. Then run ADMIXTURE in projection mode ("cp admixture\_data.5.P my\_data.5.P.in; admixture -P my\_data.bed 5"). my\_data.5.Q will then have the ethnicity estimates you're looking for.
> 8 votes

---
Tags: sam, vcf, genomics, human
---