Is there any tool for fuzzy sequence matching? Note: this question was also asked on biostars

are you aware of any tool that is able to perform error-tolerant pattern-matching search on protein FASTA files?

For example, I want to know, which proteins in my fasta file to matchADNG..C.G regexp (which represents ADNGCG pattern). However, I want to be tolerant to matching errors, meaning that I'm good with any protein that differs in any 1 letter in the motif: ADNG..C.D, ADNG..D.G, MDNG..C.G, etc are all good. Running all possible variants through grep is possible, but exponentially long for longer patterns (I usually have 15-20 letters in pattern, and scan ~10e6-7 sequences).

I am aware of the tool agrep (docs. As far as I know, it's not supported anymore, and also does not distinguish letters in which it's allowed to make error. Also, it does not support long enough patterns (with more than 9 errors -- and yes, I tried to recompile it with myself from here and could not get it work). Also, it's not designed for proteins specifically.

I am also aware of ScanProsite and Protein Pattern Find services. The first one is not error-prone, the second one is seems to be a web interface for grep.

Python solution:
e<=2 will allow up to 2 insertions, 2 deletions, or 2 SNPs in any 2 combinations

import regex
seq_list=["GGGAGGG", "GGGAAGGG", "GGGAAAGGG", "TTTATTT", "TTTAATTT", "TTTAAATTT", "GGGG", "GGG"]
for sequence in seq_list:
    if regex.search("(GGGGGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAGGG
GGGAAGGG
TTTATTT
TTTAATTT
GGGG

Can also include a wildcard

seq_list=["GGGAAAAAGGG", "GGGAAAAATGG", "GGGAAAAATTG", "GGGAAAAATTT"]
for sequence in seq_list:
    if regex.search("(GGG.*GGG|TTTTTT){e<=2}", sequence):
        print(sequence)

GGGAAAAAGGG
GGGAAAAATGG
GGGAAAAATTG


pysam or piping samtools view to a python script Is it faster to use the PySam package to run a python script on a bam

for read in samfile.fetch('chr1', 100, 120):
 print read

compared to using a pipe and reading line by line from stdin in a python script.

samtools view example.bam | python script.py.

for line in tqdm(sys.stdin, total=number_of_lines)

It's a bit hard to say with certainty, though I would suspect that offloading the BAM decompression by using a pipe will be very slightly faster. Note that decompressing and parsing the BAM file will not be the bottleneck in your processing, rather the python script itself will be.


Script to allow gene set enrichment analysis of 10x genomics data in R I have 10x single cell RNA seq data. Which R package is best suited for analysis of the 10x data matrix. What is the script to prepare the data for downstream GSEA analysis.

I have already processed samples for single cell data (for multiple other purposes). my intention is to focus on subpopulations i.e. cell type x at different time points i.e. time 0, 1, 2 etc and compare the enriched gene pathways in this one cell type at different time points. there are three files that are used (at least in the R Seurat package) with the following file extensions: .mtx, .tsv (barcodes.tsv, genes.tsv and matrix.mtx). the R Seurat package merges these files to create a "Seurat object" which helps mainly with identifying individual diffl gene expression in each barcode cell

I found a package for R: goseq that at least has been used for bulkRNA seq but I was hoping someone has used this package specifically for single cell RNA seq data, and what script has been used to merge the 3 files to create an object for analysis. pardon my terminology, but I am a novice to R and have experience primarily through package use.

Seruat will give you a list of genes which it thinks are upregulated in a particular cluster. Look at the functions that talk about marker genes - these functions basically do a DE analysis of the genes in one cluster compared to the others.

Then take that list and feed it to any standard GO analysis tool. Have a look at the topGO topKEGG and geneSetTest functions in the limma package, the GOStats package and the gsea. All should be suitable.  The GOSeq package is designed to compensate for gene length bias in RNA-seq experiments. As 10X only samples 3'-tags, there shouldn't be any gene-length bias in the data, so this shouldn't be an issue.


Confusion related to effective population size I am making a demographic model for African American population in which I say 80% of ancestry comes from African population and 20% comes from European population. I am using the time/population time from Gazave et al. (2014).

The authors doe not mention the effective population size (Ne) of African American population after admixture. Is there a publication someone can point me to that mentions the effective population size of African American (AA) or even European American (EA). I know there is not much difference between the two. It's just I need a figure to have my simulation up and running.

Note Ne is an estimate of the people size based genetic diversity and migration and is associated with the Fst statistics

That is a fantastic question and in my opinion a really important question in context. The Ne of Africa will be much bigger than the European population, however I'm less sure about the West African population which gave rise to the Afro-American population. I suspect it is still large because of the "out of Africa" theory.

My advice is to start by taking a complete look at the entire subject, from Tenesa et al (2007).  Ne is a really important parameter in any population study and extremely useful. I have feeling this could be complex, but I cannot hypotheses any further because I don't do human genetics (except for immunology).


How to change sequence format in my alignment file? I have fasta file with alingned several sequences (from MUSCLE), when I open it (e.g. notepad++) they look like:

And I want dot format of identities like  and then save it in txt file.

I failed using MEGA, BioEdit, Jalview and CodonCode, how can I do this?

This is a comment rather than answer - but not enough rep

I don't think you'll find any tool to do this exactly as you want  - how is it going to decide what's a dot and what's a letter (or a space,N or whatever) if you have multiple ambiguities or even an even number of alignments split 50:50?

Wouldn't be hard to code something which prints the whole column if it consists of non-identical nucleotides, else dots. Would be easy in R, could probably even do it in Excel without too much effort and no need to use VBA.


Sequencing center returns data in several files We send some samples to sequence and we got several (fastq.gz) files for each sample. The files are distributed at two or three folders with different dates (more than a week apart). The dates of the folders where the data is in:

2019-08-12
2019-08-21
2019-08-21
2019-09-04
2019-09-14

I asked the facility about why there are the same file in different folders, and they told us that:

  we have split sequencing of that project on different flowcells, so all
  files are correct and you have to merge the individual (duplicate) files
  to get the requested total read numbers.

How can different flow cells be from the same experiment and be run with more than a week of difference? (I thought that all the flow cells are run at the same time in each sequencing)

Can be there some batch effect? How can I check if there is batch effect if we merge the files before trimming, mapping and counting?

Usually, a single sequencing library is created from each sample. Assuming you are using barcoded indexes to identify multiple samples, the index is attached to each library, and the libraries are pooled to make a library pool.  A library pool can be frozen, and is very stable.

A pool can be sequenced multiple times by pipetting a portion of the pool and supplying it to the sequencing process.  Each time, a sample of the molecules in the pool will be sequenced.

Generally, it is true that sequencing runs are very repeatable, and that you will not suffer any batch effects.

If you want to check for batch effects you would need to analyse all of your sequence files separately, and check this first, or in parallel, to your main analysis.

One thing you might want to check is whether the same library pool was sequenced multiple times, or whether separate pools or libraries were made for each sequencing run.  In the second case, there may be more likely to be batch effects.


Continuing analysis based on a log-likelihood ratio? (Comparing methylation between two PoolSeq phenotypes) I'm trying to get some pointers on where to start on the analysis of my dataset. First, some basic info:
Background:
We performed Nanopore sequencing of two poolseq samples called UU and DD referring to damage vs undamaged. We were able to obtain approximately 3x coverage (this is currently more of a proof of concept than a hardcore association job, so the coverage is unfortunate but worth pursuing). Between these samples, we know that one is going to have different methylation patterns in response to a phenotypic trait response, but we don't know what regions along the genome are going to be differently methylated. The goal is to compare what's methylated (as in higher or lower levels of methylation) between the two samples, knowing there ought to be a lot that cancels out and come up with ideally a short list of significantly different regions we can associate to the plastic response.

Output of initial methylation calling data:
Using the raw signal level data against the reference genome and sorted, mapped BAM files the Nanopolish software was used to call methylation. This produces a dataset with the following headers:
1. Chromosome, 2. +/- Strand 3. Start vs. 4. End of the CG dinucleotide (the difference between these is 1 or two bases) 5. The read ID the dinucleotide came from 6. the log likelihood ratio calculated from the model embedded in the methylation calling process using 7. log likelihood that it's methylated and 8. log likelihood that it's unmethylated. 8. number of calling strands (always appears to be 1) 9. number of motifs (don't know about this) and finally 10. the sequence surrounding the dinucleotide, which is always short but variable in length.

I can follow the basic interpretation of the log-likelihood... if it's a positive number, there is statistical support that it's methylated. This would probably be about as much as we need, short of a p-value, if we were looking for a specific region, but we are trying to analyze the differences in de/methylation across the entire methylome. How would you suggest using the statistic values provided in our dataset to compare and visualize the significant different modifications of the dinucleotides between strands/samples?

I assume you don't want to spend 1-2 years developing your own method based on nanopore data, so the path forward is to convert the output into a format used by standard DMR calling programs. There's no standard format for that, annoyingly, but the general things you'll need are:

Coordinates of a C/CpG with coverage
Number of reads supporting methylation at that site.
Number of reads not supporting methylation at that site.

Typically you would do some filtering before/after making such a table so you only have decent quality calls (that would translate into some reasonable log likelihood for your data) and possibly only sites with decent (5x, 10x, 30x, etc.) coverage. You'd then feed that for all of the replicates in all of your groups into programs like bsseq, DSS or Metilene (to not so randomly select 3 examples).

If there aren't DMR finding programs already developed specifically for nanopore data I expect someone is developing it. The technology has a higher error rate, but given that bisulfite conversion as required in Illumina data also leads to a higher error rate that's not necessarily a huge issue if the methylation calling itself is reasonably accurate and the depth is high enough. Note that you probably need 10X coverage as a bare minimum unless your effect size is very very large.


Can I use BioEdit to quantify the sequence similarity? I am using BioEdit. It can perform an align of multiple sequences. Can this sequence similarity be outputted as a value, e.g. 70%?

No Bioedit is an alignment editor only. It is not in anyway an analysis program.

Dump the format in fasta or mega (its the option is available) and import it into Mega X here. You can quantify the diversity as a pairwise matrix, population analysis (if you have a clear a priori that you can subdivide your sequences into different groups) or phylogeny ... which provides diversity as mutations per alignment position. Its pointy and clicky and easy to use.


How to run HYPHY on multiple files non-interactively? I would like to test for positive selection in a large set of genes.

I want to have a yes/no answer to the question if gene was evolving under positive selection. Therefore I chose model BUSTED for my analysis, which is implemented in HYPHY.

When I launch HYPHY in the command-line, I get a series of question I  have to answer in the interactive form.

Is there way to perform a batch analysis involving multiple datasets, i.e. many sequence alignments?

TL;DR

I made the following bash script to generate an input file for HYPHY.

#!/bin/bash
cat << EOF
inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="$(readlink -f $1)"; // codon data
inputRedirect["03"]="$(readlink -f $2)"; // tree
inputRedirect["04"]="${3:-All}"; // Test for selection on a branch
inputRedirect["05"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);
EOF

If you want to test all the branches jointly run:

./gen_busted.sh alignment.phy tree.nwk > script.bf

Or for a particular branch:

./gen_busted.sh alignment.phy tree.nwk branchLabel > script.bf

And then run:

HYPHYMP script.bf

Now you can use the script to generate a .bf file for every sequence alignment you have, and run HYPHY for every file.

Longer version

I will describe how to make a similar script for any model, such as MEME, RELAX or PARRIS.

First, find the .bf script which performs the analysis. On GNU/Linux it should be at HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles (it might be different between different OS). In case of BUSTED your script has a name BUSTED.bf.

Now run this model in the interactive mode in order to record all the inputs required by the model. Make sure to specify full paths to all the files.

HYPHYMP HYPHY_INSTALLATION_PATH/lib/hyphy/TemplateBatchFiles/BUSTED.bf

In case of BUSTED the inputs are:

Genetic code.
Codon data (sequence alignment in phylip format).
Tree (in newick format).
Which branches to test.

Now for every input alignment you need to generate a batch file which looks like this:

inputRedirect = {};
inputRedirect["01"]="Universal"; // genetic code
inputRedirect["02"]="/path/to/alignment.phy"; // codon data
inputRedirect["03"]="/path/to/tree.nwk"; // tree
inputRedirect["04"]="All"; // Test for selection on all branches
inputRedirect["05"]="BRANCH1"; // Test for selection on branch1
inputRedirect["06"]="BRANCH2"; // Test for selection on branch2
inputRedirect["07"]=""; // complete selection

ExecuteAFile (HYPHY_LIB_DIRECTORY+"TemplateBatchFiles/BUSTED.bf", inputRedirect);

Now you can create a script which generates the file. The most important thing, don't forget to use full paths in the .bf file.


How to filter info score post-imputation? this is a very basic question, but I cannot find it explicitly stated anywhere: When exactly should I filter the info scores after imputation?

I did imputation in Impute2 and am planning on using SNPtest to analyze, but cannot tell if I should filter the scores in Impute2, SNPtest, Gtool, QCtool, or through python.

I'm hoping to get rid of the poorly imputed results (with an info score cutoff <0.5) in a population of ~4000 (half cases, half controls). After this, I'm running an association analysis -- probably through PLINK, but I'm still figuring out the best way to go from Impute2 to PLINK.

Any help would be greatly appreciated!!

In simple GWAS setups, each SNP is analyzed independently.  In those cases you can filter out SNPs with poor INFO scores at any point.

For analyses that combine information across SNPs (for example FaSTLMM), I would recommend filtering out SNPs before running the step that combines information across SNPs.

So, if you filter right after running Impute2, you should be safe.


Check if a string is ensembl gene id or not in python I have many scRNAseq datasets saved as hdf5 files (hierarchical data format*), and when I open them keys sometimes are confused, so I need a way of determining from string values themselves whether the array contains gene symbols or ensembl ids. What would be a regular expression in python to check for it?
*, A database format for quickly accessing large amounts of data, which is accessed using keys. The format is quicker than pickle/unpickle and can be viewed via a hdf5 viewer, a bit of info is here

Fixed it myself actually, not that hard. I found regular expression for ensembl ids here:

https://www.biostars.org/p/99142/

And used the re package in python:

import re

ensembl_regexp = 'ENS[A-Z]+[0-9]{11}'

is_ensembl = np.all([bool(re.match(ensembl_regexp, gene)) for gene in gene_symbols])
    if is_ensembl:
        gene_symbols = list(f['ensembl_ids'])

Where gene_symbols is a numpy array that I need to check whether it has ensembl ids, and if yes, switch to using gene symbols which are going to be under ensembl_ids key due to a mistake of someone who made the hdf5 dataset.


dbnSNP frequency anomalies Sometimes dbSNP reports very different allele frequencies for different large-scale genome projects e.g. between 1000 Genomes and GnomAD
rs11822440 1000Genomes A=0.4629    C=0.5371 GnomAD A=0.99997   T=0.00003
rs28531832 1000Genomes T=0.4379    C=0.5621 GnomAD T=0.99993   A=0.00007
rs78931618 1000Genomes A=0.7334    C=0.2666 GnomAD A=0.9999    T=0.0001
rs75716913 1000Genomes T=0.9006    C=0.0994 GnomAD T=0.99149   A=0.00851
Is this an issue with one or other of these projects data processing? Is something else going on here biologically or technically that might cause this?

I'll stick my neck out and simply state that allopatry alone would account for the pattern you are observing. Thus populations can exihibit variation in allelic profiles, particularly where Mendellian inheritance is operating (i.e. higher eukaryotes), singly due to geographic isolation.
This represents the Meyers view of population genetics and would account for variation within neutral theory, e.g. silent mutations. Obviously I'm guessing the populations were selected by geographic foci in generating the allelic frequencies.


Is there a standard method to validate whether a bigWig file is correct? I have a bigWig file of extension *.bw (which is standard), but I am suspicious that this file may be of an incorrect format.

Is there a standard method in the field to validate a *.bigWig file?

There's no standard validation method, but a simple one would be python -c "import pyBigWig; bw = pyBigWig.open('your_file.bw'); print(bw.IsBigWig());". If you get True then it's a bigWig file. If you get False or an error then it's not (False indicates a bigBed file).

Alternatively, file your_file.bw should return data.


Why SeqIO.parse method isn't working? I am trying to follow the tutorial from from the Biopython website here  and I am right at the beginning (2.4.1), where I am trying a Simple FASTA parsing example which is
from Bio import SeqIO
for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
    print(seq_record.id)
    print(repr(seq_record.seq))
    print(len(seq_record))

For an unknown reason, I get the following error:
---------------------------------------------------------------------------
FileNotFoundError                         Traceback (most recent call last)
<ipython-input-19-c64128d7dcfd> in <module>
      1 from Bio import SeqIO
----> 2 for seq_record in SeqIO.parse("ls_orchid.fasta", "fasta"):
      3     print(seq_record.id)
      4     print(repr(seq_record.seq))
      5     print(len(seq_record))

~\Anaconda3\lib\site-packages\Bio\SeqIO\__init__.py in parse(handle, format, alphabet)
    605     iterator_generator = _FormatToIterator.get(format)
    606     if iterator_generator:
--> 607         return iterator_generator(handle)
    608     if format in AlignIO._FormatToIterator:
    609         # Use Bio.AlignIO to read in the alignments

~\Anaconda3\lib\site-packages\Bio\SeqIO\FastaIO.py in __init__(self, source, alphabet, title2ids)
    181             raise ValueError("The alphabet argument is no longer supported")
    182         self.title2ids = title2ids
--> 183         super().__init__(source, mode="t", fmt="Fasta")
    184
    185     def parse(self, handle):

~\Anaconda3\lib\site-packages\Bio\SeqIO\Interfaces.py in __init__(self, source, alphabet, mode, fmt)
     45             raise ValueError("The alphabet argument is no longer supported")
     46         try:
---> 47             self.stream = open(source, "r" + mode)
     48             self.should_close_stream = True
     49         except TypeError:  # not a path, assume we received a stream

FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'

I am using a Jupyter Notebook.

Personally I think its late and we're all tired.

"FileNotFoundError: [Errno 2] No such file or directory: 'ls_orchid.fasta'"

The traceback is pointing to line 2, the SeqIO as you point out.
In my opinion, the script should be in the same directory as the fasta file. You know how to identify the path in bash (pwd or echo $PWD), although you must be using Windows (Anaconda has \ rather than /). Alternatively the file names has a typo. You know the
for seq_record in SeqIO.parse("~\path\ls_orchid.fasta", "fasta"): # where path is the dir(s) leading to ls_orchid.fasta, but obviously use / if its Linux.

I would assume you can alternatively dump the .ipynb file in the location where the fasta file is.
As personal choice I prefer a formal IDE such as 'Visual Studio Code' or "PyCharm". In my personal opinion, I am aware some will disagree, Jupyter notebook good is for running code on a remote commercial server for data science work. However, the traceback - which you clearly have activated - looks good. I can see that you're Anaconda3 installation looks good, you must be using Windows (which I've never worked with) and you're running Jupiter from Anaconda3 ... so cool.


efficient counting of dinucleotides/trinucleotides on fastq reads? What's an efficient way of counting of dinucleotides/trinucleotide pattern on fastq.gz file reads?

I know there are tools like seqtk that will be very efficient at reading through the .fastq.gz files, and tools like cutadapt that are good at finding long(ish) adapter patterns in the fastq.gz PE of files, but, what would be the best tool to count the number of reads where a di-/trinucleotide is present N-number of times in the read?

E.g.

countingtool --pattern CGC file_R1.fastq.gz

Result: column two is number of times the pattern is seen (exclusive numbers, no overlap)

CGC    0     51341234
CGC    1     13948130
CGC    2      1234344
CGC    3        51344
CGC    4         1343

Well, I'm sure you can find a more efficient approach, but here's a simplistic solution using standard UNIX tools:

zcat file.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

I tested this on a 2.6G fastq file (exome data) and it took:

real    1m18.987s
user    1m41.392s
sys     0m10.684s

On a 66G WGS sequence file, it took:

real    34m6.771s
user    40m16.452s
sys     4m10.176s

The extra grep CGC isn't needed, you can do the match in awk, but it does speed things up a little bit. You can always get the reads with 0 matches by counting the number of reads and subtracting the ones with 1. Alternatively, remove the grep:

zcat file.fq.gz | awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
    sort -rnk2

Following the OP's suggestion, I also tested with pigz and, it turns out, pigz can fly!

time unpigz -c exome.fastq.gz | grep CGC |
    awk -F'CGC' 'NR%4==2{a[NF-1]++}END{for(i in a){print i,a[i]}}' |
         sort -rnk2

real    0m44.012s
user    1m19.544s
sys     0m14.900s

And for the 66G WGS:

real    17m29.935s
user    29m23.816s
sys     5m17.856s


Uniprot: Determining protein and organism given a DNA sequence I would like to find the identity of a protein, together with the organism name and kingdom.

The given DNA sequence is below:

ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTATGATATTCTTC
TTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATTCCAGATGAGTGTACCGGAGT
TAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTTAGTATGTGGTCTACTTTTGATTACGGGTTC
ATCTTGGATGGAGGAAGGTTCAGGAACGGCCTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCAT
AGCGGACTTGCTGTAGATACGTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCA
TCAACCTTATATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC
CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGCTATCACGATG
CTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGGCGGTGATCCTATGTTGTATC
AGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATATCATCATACTTCCAGTATTCGGTATAGTCAG

I've been asked to translate this DNA sequence into the 6 possible frames which are:

5'3' Frame 1
TVVR-SSVVQRPNNITCNAYDILLHHARHDGRTRQLTSAIPDECTGVSIPKD--HRYMIF
SMWSTFDYGFILDGGRFRNGLNRLSTTSAHCKS-RTCCRYVHYRIAHGRCKLPYRKHQPY
MYNRLCPPFTHGDAAVITLSLIHYNHCSVTHRSCACASRCYHDATH--KLEYQLL-QFGR
R-SYVVSALILGVWASRSLYHHTSSIRYSQ
5'3' Frame 2
LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL
VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV
5'3' Frame 3
CCSVIISCTTS-QHHMQCL-YSSSSCQARWQD-ATY-CHSR-VYRS-YSQRLITSVYDF-
YVVYF-LRVHLGWRKVQERPEPSIHH-RSLQVIADLL-IRSLSHCTWPVQAPLQEASTLY
VQSPMPAVHSWRCCSHHFIPDPLQSLQRYS-ELCLC-QVLSRCYSLIEVGVPASMTVRQA
VILCCISTYSGCLGIQKSISSYFQYSV-S
3'5' Frame 1
LTIPNTGSMMI-TSGCPNTQNKC-YNIGSPPAELS-KLVLQLLSVSSIVIAPASTGTTPM
SNAAVIVMDQG-SDDCSIAMSERRA-AIVHIRLMLPVRELAPAMCNAIMNVSTASPL-LA
VSASGG-TVQAVPEPSSIQDEPVIKSRPHTKNHIPMLLIFGNTNSGTLIWNGTSKLPSPA
IVPGMMKKNIISIACDVVRTLYN--SPNNS
3'5' Frame 2
-LYRILEV--YRLLDAQTPRISADTT-DHRLPNCHRSWYSNFYQ-VAS--HLLAQAQLL-
VTLQ-L-WIRDKVMTAASP-VNGGHRRLYI-G-CFL-GSLHRPCAMR--TYLQQVRYDLQ
-ALVVDRRFRPFLNLPPSKMNP-SKVDHILKIIYRCY-SLGILTPVHSSGMALVSCLVLP
SCLA--RRIS-ALHVMLLGRCTTDDHRTT
3'5' Frame 3
DYTEYWKYDDIDFWMPKHPE-VLIQHRITACRTVIEAGTPTSISE-HRDSTC-HRHNSYE
-RCSDCNGSGIK--LQHRHE-TAGIGDCTYKVDASCKGACTGHVQCDNERIYSKSAMTCS
ER-WWIDGSGRS-TFLHPR-TRNQK-TTY-KSYTDVINLWEY-LRYTHLEWH--VA-SCH
RAWHDEEEYHKHCM-CC-DVVQLMITEQQ

I ran jobs on Uniprot for all these sequences. Here are the results:

www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E0FZ
www.uniprot.org/blast/uniprot/B2017100383C3DD8CE55183C76102DC5D3A26728BEBF5BAE
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082E5FP
www.uniprot.org/blast/uniprot/B20171003A7434721E10EE6586998A056CCD0537E082F08D
nothing
www.uniprot.org/blast/uniprot/B20171003F725F458AC8690F874DD868E4ED79B88CE32C4B

But how can I decide the identity of the protein and the name of the organism?

In the first one, the organism is Pacific ridley turtle  but in the second one it is Diplonema sp. ATCC 50224. So I can only decide at a kingdom level, but not to the organism, right?

Well, let's see. The first step is to simply run a nucleotide blast (blastn) against the nr (non-redundant) database. I did this using the NCBI blast tool. While, as expected, I got many hits, there was only one hit which matched the entire length of the query sequence perfectly:

Diplonema sp. ATCC 50224 cytochrome c oxidase subunit 1 (cox1) mRNA, complete cds; mitochondrial
Sequence ID: JF698657.1
Length: 1612Number of Matches: 1
Range 1: 141 to 770
Alignment statistics for match #1
Score   Expect  Identities  Gaps    Strand
1164 bits(630)  0.0 630/630(100%)   0/630(0%)   Plus/Plus
Query  1    ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  60
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  141  ACTGTTGTTCGGTGATCATCAGTTGTACAACGTCCTAACAACATCACATGCAATGCTTAT  200

Query  61   GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  120
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  201  GATATTCTTCTTCATCATGCCAGGCACGATGGCAGGACTAGGCAACTTACTAGTGCCATT  260

Query  121  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  180
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  261  CCAGATGAGTGTACCGGAGTTAGTATTCCCAAAGATTAATAACATCGGTATATGATTTTT  320

Query  181  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  240
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  321  AGTATGTGGTCTACTTTTGATTACGGGTTCATCTTGGATGGAGGAAGGTTCAGGAACGGC  380

Query  241  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  300
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  381  CTGAACCGTCTATCCACCACTAGCGCTCACTGCAAGTCATAGCGGACTTGCTGTAGATAC  440

Query  301  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  360
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  441  GTTCATTATCGCATTGCACATGGCCGGTGCAAGCTCCCTTACAGGAAGCATCAACCTTAT  500

Query  361  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  420
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  501  ATGTACAATCGCCTATGCCCGCCGTTCACTCATGGCGATGCTGCAGTCATCACTTTATCC  560

Query  421  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  480
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  561  CTGATCCATTACAATCACTGCAGCGTTACTCATAGGAGTTGTGCCTGTGCTAGCAGGTGC  620

Query  481  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  540
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  621  TATCACGATGCTACTCACTGATAGAAGTTGGAGTACCAGCTTCTATGACAGTTCGGCAGG  680

Query  541  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  600
            ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Sbjct  681  CGGTGATCCTATGTTGTATCAGCACTTATTCTGGGTGTTTGGGCATCCAGAAGTCTATAT  740

Query  601  CATCATACTTCCAGTATTCGGTATAGTCAG  630
            ||||||||||||||||||||||||||||||
Sbjct  741  CATCATACTTCCAGTATTCGGTATAGTCAG  770

So, we now know the sequence is from the gene for mitochondrial Cytochrome C Oxidase subunit 1 from Diplonema sp. ATCC 50224. That the sequence is mitochondrial, is very important since the mitochondrion of this species has a slightly different genetic code. If you go to the NCBI's page for this sequence (JF698657.1), you will see this line:

/transl_table=4

That tells you it is using translation table 4:

4. The Mold, Protozoan, and Coelenterate Mitochondrial Code and the Mycoplasma/Spiroplasma Code (transl_table=4)

    AAs  = FFLLSSSSYY**CCWWLLLLPPPPHHQQRRRRIIIMTTTTNNKKSSRRVVVVAAAADDEEGGGG
  Starts = --MM------**-------M------------MMMM---------------M------------
  Base1  = TTTTTTTTTTTTTTTTCCCCCCCCCCCCCCCCAAAAAAAAAAAAAAAAGGGGGGGGGGGGGGGG
  Base2  = TTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGGTTTTCCCCAAAAGGGG
  Base3  = TCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAGTCAG

See how all of your translated reading frames have in-frame stop codons? That's because you translated them using the standard genetic code instead of the correct mitochondrial one. As you can see above, in this code, the codon TGA is not a STOP, as in the standard code, and instead codes for tyrosine (W).

So, if you look at the result of the protein blast you ran on UniProt, you will see:

Query1   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI-FL   59
         LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGI FL
F5CI61   LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL   107

Query60  VCGLLLITGSSWMEEGSGTA-TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   118
         VCGLLLITGSSWMEEGSGTA TVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
F5CI61   VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI   167

Query119 CTIAYARRSLMAMLQSSLYP-SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   177
         CTIAYARRSLMAMLQSSLYP SITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
F5CI61   CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG   227

Query178 GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  206
         GDPMLYQHLFWVFGHPEVYIIILPVFGIV
F5CI61   GDPMLYQHLFWVFGHPEVYIIILPVFGIV                                  256

The only locations where there isn't a 100% identical match is where a W is aligned to a STOP (- in the Uniprot output). So, your actual protein sequence is the 2nd frame of the 5' to 3' direction which, when translated using the correct code, is:

LLFGDHQLYNVLTTSHAMLMIFFFIMPGTMAGLGNLLVPFQMSVPELVFPKINNIGIWFL
VCGLLLITGSSWMEEGSGTAWTVYPPLALTASHSGLAVDTFIIALHMAGASSLTGSINLI
CTIAYARRSLMAMLQSSLYPWSITITAALLIGVVPVLAGAITMLLTDRSWSTSFYDSSAG
GDPMLYQHLFWVFGHPEVYIIILPVFGIV


differentially expressed intron analysis? I have used Deseq for estimating the differentially expressed introns for the RNA seq data using the UCSC intron bed file. can Deseq be used for differentially expressed intron analysis?

Typical tools to use for this would be the DEXseq packages (as done here) or rMATS, which will calculate intron retention rates.

I would generally suggest not summing intronic counts by gene and using that in DESeq2 (I assume that's what you meant by DESeq, since there's no reason to use the original DESeq for anything). What you'd be getting if you did this is less expressed introns and more sample-specific bias in enriching for mature versus pre-mRNA.


How to apply upperquartile normalization on RSEM expected counts? I see that TCGA RNASeq V2 RSEM data is normalized with upper-quartile normalization.

After doing Quantification with RSEM with the samples I have, I got "genes.results" as output which has gene id, transcript id(s), length, expected count, and FPKM. So, from all the sample output files I got the gene_id and expected_counts [from all samples]. for eg: it looks like following:

Gene_id S1  S2  S3  S4  S5
ENSG00000000003 1355    1121    242 709 1071
ENSG00000000005 5   0   0   0   0
ENSG00000000419 839 1345    48  976 1234
ENSG00000000457 429 1803.08 386 1156    628
ENSG00000000460 404 1444    523 479 1260
ENSG00000000938 294 312 93  189 683
ENSG00000000971 3911    4633.92 264 2863    5194
ENSG00000001036 1546    2276    271 1641    2141
ENSG00000001084 1945    2125    490 980 1533
ENSG00000001167 2054    4069.72 3457    2075    2207
ENSG00000001460 153 339 77  221 273

I want to apply "upper-quartile" normalization on this data. But doesn't know which R package I can use and which function?

Is anyone aware about this?

P.S. This is not for differential analysis

You can use the quantile function in base R to get the value of a particular quantile (e.g. 0.75 for the upper quartile). This can then be used as a factor for normalisation: divide the observed expression by this number.

## set up simulated data
data.mat <- floor(10^data.frame(sapply(1:5,function(x){rnorm(100)})));
colnames(data.mat) <- paste0("S",1:5);
data.mat[data.mat > 0] <- floor(data.mat[data.mat > 0]);

It's generally a good idea to first remove genes that have no counts in any samples. Rows can be removed using a logical selection based on the total expression in any given row.

## Only retain rows that have a positive total expression per gene
data.mat <- data.mat[rowSums(data.mat) > 0,];

Two methods of upper quartile determination are presented here, one which is based on the total counts, and one based on the expressed counts (i.e. excluding values with no expression). These both use the quantile function of R. The apply function calculates statistics that are based on a single specified dimension (or multiple dimensions) of the data. The first dimension of a 2D matrix in R corresponds to the row, while the second corresponds to the column. By using 2 in the apply function, data are collected along the second dimension (i.e. the column), and statistics generated from that:

data.quantileAll <- apply(data.mat, 2, function(x){quantile(x, 0.75)});
data.quantileExpressed <-
  apply(data.mat, 2, function(x){quantile(x[x>0], 0.75)});

Normalisation is applied by dividing the values by the column-based upper quartile statistic. R will repeat values along the last dimension first (which in this case is the column), so the matrix needs to be transposed first (with columns as genes) to get the arithmetic working properly, with a final transposition to return the orientation back to the original direction (with columns as samples).

data.norm <- t(t(data.mat) / data.quantileExpressed);


low-memory high-speed replacement for Picard MarkDuplicates I am running Picard MarkDuplicates with the following parameters below. On the file described, it takes about 41.6Gb of RAM memory and about 20-25 minutes to compute (only uses 1 core AFAICS).

java -Xmx54890m -jar picard.jar MarkDuplicates \
    I=/home/dnanexus/in/sorted_bam/TST73-86-2IC_S16_L00.bml.GRCh38.karyo.bam \
    O=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.deduplicated.bam \
    M=./TST73-86-2IC_S16_L00.bml.GRCh38.karyo.duplication_metrics \
    CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT REMOVE_DUPLICATES=true

What equivalent free/open source software (samtools?, sambamba?) can I use to obtain the same deduplicated bam file but at a lower cost, considering the instance used here has 8 cores but only one of them is being used, and has about 68Gb of RAM out of which it's using 46Gb odd at peak.

I think previous versions of samtools/sambamba didn't behave the same way as Picard MarkDuplicates with the parameters above, but I would like to know which recent versions would give me the same output.

As you suggested, sambamba is faster at marking duplicates than picard (it's also multithreaded). Recent versions of samtools have a rewritten duplicate marking algorithm, though I doubt it'll be as quick as sambamba. 46GB of RAM seems excessive for marking duplicates unless you're having it store the whole file in memory.


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

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

Any ideas?

Thank you!

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

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


A tool or webserver for building PSSM matrix I have some protein sequences and I want to build a position-specific scoring matrix (PSSM) for them and then upload this PSSM to NCBI PSI-BLAST. I used CHAPS program for this pupose but uploading the output PSSM gave me an error in NCBI PSI-BLAST. Do you know any tool or webserver for getting PSSMs for a group of protein sequences which then can work in NCBI PSI-BLAST?
This is a part of my PSSM file:

PssmWithParameters ::= {
pssm {
    isProtein TRUE ,
    numRows 28 ,
    numColumns 131 ,
    byRow FALSE ,
    query
    seq {
        id {
            other {
                accession "WP_000208753" } } ,
        inst {
            repr raw ,
            mol aa ,
            length 131 ,
            seq-data
            ncbieaa "MTTKRKPYVRPMTSTWWKKLPFYRFYMLREGTAVPAVWFSIELIFGLFALKNGPEAW
AGFIDFLQNPVIVIINLITLAAALLHTKTWFELAPKAANIIVKDEKMGPEPIIKSLWAVTVVATIVILFVALYW" } } ,
    intermediateData {
        freqRatios {
            { 0, 10, 0 } ,
            { 564418841, 10, -10 } ,
            { 0, 10, 0 } ,
            { 11768571, 10, -9 } ,
            { 185838265, 10, -10 } ,
            { 31496547, 10, -9 } ,
            { 3872857, 10, -8 } ,
            { 291750464, 10, -10 } ,
            { 128450763, 10, -10 } ,
            { 759856221, 10, -10 } ,
            { 359173937, 10, -10 } ,
            { 179865517, 10, -9 } ,
            { 14537895, 10, -8 } ,
            { 212921456, 10, -10 } ,
            { 220554141, 10, -10 } ,
            { 368516324, 10, -10 } ,
            { 319343525, 10, -10 } ,
            { 426173953, 10, -10 } ,
            { 463659523, 10, -10 } ,
            { 817322186, 10, -10 } ,
            { 811693964, 10, -11 } ,
            { 0, 10, 0 } ,
            { 227810064, 10, -10 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 0, 10, 0 } ,
            { 768325726, 10, -10 } ,
            { 0, 10, 0 } ,
            { 1425562, 10, -8 } ,
            { 372685285, 10, -10 } ,
            { 466727421, 10, -10 } ,
            { 185741084, 10, -10 } ,
            { 427328646, 10, -10 } ,
            { 122594965, 10, -10 } ,

If you take a look at my answer in this BioStars post, you can generate a PSSM using AlignIO in Biopython:


Putting cells in Seurat dimplot selectively I am going to adjust Seurat dimplot in a way avoiding some cells so both my dimplot and heatmap look nice

If you please consider this picture, you would see some cells are far from the clusters so I want to avoid them in dimplot and of course for heatmap (coming from finding markers).

[![enter image description here][1]][1]

Can you suggest me a piece of code doing that?
Thanks

If you're using a GUI you could select the cells interactively:

plot <- DimPlot(seurat_obj, reduction = "umap")

Then select the cells by clicking around them

select.cells <- CellSelector(plot = plot)
Idents(seurat_obj, cells = select.cells) <- "SubCells"

and subset based on these cells

sub_cells <- WhichCells(seurat_obj, idents = "SubCells")
sub_obj <- subset(seurat_obj, cells = sub_cells)

You could then use subobject in your downstream dimplot.

See here for some of the interactive options. I guess you could outright subset on UMAP coordinates, but that may not always be easy if your boundary isn't linear.


What is the best way to map/align my reads on a given genome? I am frequently using a ballgown package for my rnaseq analysis, but recently I have had a new task to have my reads mapped on two different genomes to understand the level of alignment between the two, therefore I wonder what is the best tool to do this? I have single-end reads.
What I need is to have a percentage of uniquely mapped fragments, mapped to multiple loci and, say, unmapped. Can bowtie do the job? Is there a standard explanation procedure explained elsewhere? Thanks.

This is a pretty straightforward analysis, assuming that you want to measure something like "how many reads align to genome 1, and how many reads align to genome 2?".

The bowtie2 manual has a section on this.
There are other tutorials for bowtie2.
bwa is also well-suited to this task.

You should be able to then use samtools to compute statistics on your alignments.


Subsetting from seurat object based on orig.ident? I am pretty new to Seurat. I want to subset from my original seurat object (BC3) meta.data based on orig.ident. however, when i use subset(), it returns with Error.
ER_HER_P <- subset(BC3, idents = c("BC03"))
Error: No cells found

table(Idents(BC3))
BC01   BC02   BC03 BC03LN   BC04   BC05   BC06   BC07 BC07LN   BC08   BC09   BC10   BC11
21     50     33     51     42     74     15     48     52     21     55     13     11

There are 33 cells under the identity. I checked the active.ident to make sure the identity has not shifted to any other column, but still I am getting the error?
head(BC3@active.ident)
[1] BC01 BC01 BC01 BC01 BC01 BC01
Levels: BC01 BC02 BC03 BC03LN BC04 BC05 BC06 BC07 BC07LN BC08 BC09 BC10 BC11

If you are going to use idents like that, make sure that you have told the software what your default ident category is.
This works for me, with the metadata column being called "group", and "endo" being one possible group there.
Idents(combined.all) <- "group"
endo_subset <- subset(combined.all, idents = c("endo"))

I think this is basically what you did, but I think this looks a little nicer


How do I plot clusters based on DNA sequence alignment? I have a group of about 1050 sequences that have been aligned using blastn. I am now trying to create a cluster plot of the sequences in order to identify similar sequences, something like this:

However, I am having trouble finding a tool to do this. My professor recommended using ClustalW, and so I have spent plenty of time searching for a tool that can take a sequence alignment and plot them such as in the linked image. However, I have not been successful in finding what I need and am beginning to think I am missing something or something is going over my head. If anyone could assist with finding the tool I need I would very much appreciate it. Thanks so much!

I think one of the reasons you struggle is that clustering DNA sequences is not a clearly defined task. In general intention of clustering is to reconstruct, or approximate the relatidness of the DNA seqeunces.
If all these sequences are homologous you can do a multiple sequence alignment (using MAFFT or Clustal) and build a phylogenetic tree (using something like IQ tree). If your prof suggested ClustalW, this is probably what they meant. Evolutionary speaking, that's the optimal method for getting relations of sequences, but they must be "alignable". Also with this many sequences, it will be computationally quite heavy.
For within species relatedness, phylogenetic reconstruction is sometimes impossible, because recombination breaks down the linkage of individual alleles. However, if all of the sequences are long enough, you can do a principal component analysis on SNPs, for example using something like plink. The plot you showed actually looks a lot like PCA (it even has the explained variability on axes).
Finally, regardless of their relationship, you can always do a dendrogram based on a similarity matrix, but for that, you need to define similarity. That will be usually something kmer-based. I am not following this field very closely lately, but there are some recent papers about kmer-based clusterings, like this one.
When looking for the paper I also run into this one, which proposes genome clustering using "Inter-nucleotide Covariance". No idea how useful/sound it is. Just found the abstract quite interesting.
Long story short, it really depends on what you are after and what are your sequences! If you can build a reliable tree, it is the most powerful way to determine the relationships among sequences.


Error in if (is.na(sra_acc$run[i])) { : argument is of length zero I am trying to write a code for importing sra files as fastq files in Rstudio (for obtaining and aligning RNA-seq reads from public repositories). I keep getting "Error in if (is.na(sra_acc$run[i])) { : argument is of length zero" right after I run this code. Any suggestions on what might be going wrong? My code is as follows:
library(SRAdb)
sqlfile <-'SRAmetadb.sqlite'
if(!file.exists('SRAmetadb.sqlite')) sqlfile <<- getSRAdbFile()
sra_con <- dbConnect(SQLite(),sqlfile)
sraInf <- getSRAinfo("SRP045534",sra_con, sraType="sra")

Error in if (is.na(sra_acc$run[i])) { : argument is of length zero

After checking with file.exists(), I would suggest you try the same command with another SRA run. It may happen that some runs are not accessible for reasons I do not know. I always have problems in accessing sra formatted files and use a workaround to download fastq.
Here is my not official workaround:
runlist<-sraConvert("SRP045534",sra_con=sra_con)
runlist[1:2,]

study submission    sample experiment        run
1 SRP045534  SRA179382 SRS685546  SRX681992 SRR1552451
2 SRP045534  SRA179382 SRS685544  SRX681990 SRR1552449

for (dloadrun in runlist$run)
{
    getSRAfile(dloadrun,sra_con=sra_con,fileType="fastq")
}

This should iteratively download all the runs associated to the project you are studying.
It is possible that you can do the same with sra type file, but I have never been able to use them with sraDB, and never investigated the issue.


Using Seurat to compare mutant vs.wt I am interested in using Seurat to compare wild type vs Mutant. I don't know how to use the package. How can I test whether mutant mice, that have deleted gene, cluster together?

Single-cell analysis to compare samples is a long a difficult process. There is very good documentation for 10x Genomics cellranger, the DropSeq Pipeline and the Seurat R package. These tools all have GitHub repositories and the authors are very responsive if you encounter issues.

Depending on the technology used to generate the data, you'll need to use either cellranger or DropSeq to process the FASTQ files. These are designed to account for the the different experimental design. Still, it is important to ensure that your samples have been processed in the same way to reduce batch effects.

This question is to open-ended to cover all of the details needed to perform the analysis but here is some of the main things to consider and the general process to perform such as analysis. You will need competence using the command-line (shell) and programming in R (or Python) to perform single-cell analysis. I recommend performing this analysis on a remote server as some steps are memory intensive.

The general overview to compare samples is:

Demultiplex the reads for different indicies. You need to do this even if you've only put one sample per lane (to remove the index from the reads). cellranger mkfastq or Illumina's bcl2fastq will do this.
Obtain a reference genome (FASTA) and gene annotation (GTF) for the species you are working with. You can prepare a reference transcriptome with cellranger mkgtf and cellranger mkref.
Run cellranger count or the DropSeq pipeline on each sample separately. These will both perform STAR (splice-aware) alignment of paired-end RNA-Seq reads and count UMIs for each cell-barcode and gene. Cells will be identified and filtered by an automatic UMI threshold (this can be changed by forcing with expected number of cells). This will return a gene-barcode matrix.
You need to aggregate results of cellranger runs for different samples with cellranger aggr. This will perform downsampling by default to normalise the number of reads between samples. It will also perform graph-based (Louvian) clustering on the combined samples and return HTML and Cloupe summary data to explore, in addition to a combined gene-barcode matrix for all samples.
Single-cell RNA-Seq experiments are subject to batch effects. The CCA method from Seurat and the MNN method from scran are both available in the respective R packages to account for this (with different approaches). Batch effect correction needs to be performed before downstream analysis to ensure comparisons are valid.

This methods require an overlap between the samples to estimate technical errors and account for them. You need to bear in mind that neither of these methods works unless:

some of the cells detected in each sample are expected to be the same cell type
biological replicates of each sample have been performed

Once you have performed this correction, then standard Seurat workflows such as those shown in the tutorials can be performed which will identify clusters. You can use these to identify clusters specific to each genotype. Bear in mind that some overlap should be expected and if there isn't one, batch effect correction will over-correct for this.


Expected allele frequency distribution of SNVs in real NGS data I have a huge amount of ~20x human WGS samples, aligned, and all SNVs that were called with GATK under standard germline parameters set.

What I need to do is to model SNVs Allele Frequency (AF) for different underlying Copy Numbers. I'd better provide a toy example. For particular genomic region X:

If X is presented by 2 copies for the particular samples, we expect AF to be super-close to 1 or to 0.5.

If X is presented by 4 copies, I expect any particular AF to be close to 0.25, 0.5, 0.75 or 1.

Of course, I can use Binomial Distribution for these purposes. However, as we know, the distribution is not exactly Binomial due to alignment/sequencing biases and the median AF for all heterozygous SNVs is more close to 0.48 but not to 0.5 as we would expect. Another thing: for high copy numbers we expect higher coverages. And GATK use several filters so I suppose that we will not see SNVs with AF like 0.125 (in case if the segment has ploidy 8) - despite the super high coverage there GATK may reject this "weird" AF.

I have read several papers that model SNVs AFs (and I agree that Beta Binomial Distribution may be quite accurate), however, I was not convinced enough that I should use the particular modelling. From your experience (in case if you do SNVs calling), which probabilistic distribution should I use? How should I estimate parameters for each of them (should I expect for CN4 AF=0.5 more frequent than AF=0.75 or vice versa, how to estimate this from data)?

UPD: For simplicity we can say that we have a lot of previously identified regions with ploidy different from CN2, and I can take these coordinates from here. So I can use more or less "supervised" learning for parameters' estimation.

I don't have enough experience to answer which probabilistic distribution should be used.

However, this questions also also asks how to estimate parameters of the distributions. If a binomial distribution is chosen, then Heng Li's paper titled "A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data" 1 is probably the definitive one. Section 2.3.1 of that paper describes an EM algorithm for estimating allele frequencies from multiple samples under the assumption of Hardy-Weinberg equilibrium for arbitrary but constant ploidy.


Real time transcript profiles Do any methods exist (or are in the process of development) for investigating transcript data without lysing the cells, i.e, destroying the sample?

Answer:
There are no methods that exist for investigating transcript data inside a cell without destroying the cell. Separating RNA from cells requires essentially destroying the cells and using a chemical cocktail to split its constituent components.

Alternative 1: One option is to investigate extracellular transcripts of living cells is to centrifuge cells in culture and isolating RNA from the growth media supernatant [1]. Cell-free DNA and cell-free RNA are hot topics in contemporary biology and bioinformatics as a potential source for disease biomarkers. [2]

Alternative 2: Another option to monitor a transcript from a living cell is to use a transgenic construct with a reporter molecule like luciferase or a fluorescent protein to produce light or become fluorescent when a promotor is activated [3, 4]. This does not directly measure transcript quantity but answers a similar question: how does transcription change in a living cell?

Sources:

[1] Théry, C., Amigorena, S., Raposo, G., & Clayton, A. (2006). Isolation and characterization of exosomes from cell culture supernatants and biological fluids. Current protocols in cell biology, 30(1), 3-22.

[2] Yu, X. M., Wu, Y. C., Liu, X., Huang, X. C., Hou, X. X., Wang, J. L., ... & Ling, Z. Q. (2016). Cell-Free RNA content in peripheral blood as potential biomarkers for detecting circulating tumor cells in non-small cell lung carcinoma. International journal of molecular sciences, 17(11), 1845.

[3] Kaskova, Z. M., Tsarkova, A. S., & Yampolsky, I. V. (2016). 1001 lights: luciferins, luciferases, their mechanisms of action and applications in chemical analysis, biology and medicine. Chemical Society Reviews, 45(21), 6048-6077.

[4] Rao, J., Dragulescu-Andrasi, A., & Yao, H. (2007). Fluorescence imaging in vivo: recent advances. Current opinion in biotechnology, 18(1), 17-25.


Clarification about how mapping to a reference genome works in RNA seq When mapping reads to a reference genome, how is it possible to tell which part of the genome the read is referencing if the genome possibly has the same repeated sequences in different areas? What I mean by this is let's say for the sake of argument that you have a read "AGCTAGCT". If the reference genome has two places where the sequence "AGCTAGCT" is present, how do we tell which one it's coming from? I'm assuming that the length of the read is chosen long enough such that this doesn't happen but perhaps I'm wrong?

You can't always know exactly where every read originated, even if you use long reads to sequence full-length transcripts (odds are there's at least one transcript that's duplicated in the genome). This is also not exactly a problem, such reads can either be discarded or reported as multimapping.


what does a sequence look like before alignment? I am confused with Sequence alignment, which is

  a way of arranging the sequences of DNA, RNA, or protein to
  identify regions of similarity that may be a consequence of
  functional, structural, or evolutionary relationships between the
  sequences.

to understand this easily, I would like to use this image (image_1)

the image above (image_1) is the primary part of nucleic acid structure (primary, secondary, tertiary, and quaternary), shown below (image_2)

is the sequence gcgc... an aligned sequence? if yes, what does the sequence look like before alignment?

Sequence alignment typically refers to the primary sequence. That is, we ignore the complex 3D structure of the molecule and its spatial arrangement in the cell. Instead, we use a symbolic representation composed of symbols like A, C, G, T, and U. Alignment is a conceptual tool that, as you pointed out, is used to find regions of similarity between sequences. Therefore, it only makes sense in the context of two or more similar sequences.

By way of example, consider a simple case of two sequences X and Y. This is what they would look like before alignment. The similarities between the two sequences are unclear.

X: ATTGACATAGGCACTTCGACTA
Y: GACATCTTCGGCTAA

This is what they would look like after alignment. The similarities (and differences) are now much easier to see.

X: ATTGACATAGGCACTTCCACTA-
      |||||     ||||  |||
Y: ---GACAT-----CTTCGGCTAA

This arrangement identifies regions of similarity that may be related through evolution. This doesn't always mean perfect agreement though. As this example shows, sometimes mismatches or gaps are allowed so that long stretches of identical matches are lined up correctly. Sometimes, for very closely related sequences, the optimal alignment doesn't require any gaps at all. This seems to be the case with some of your recent questions about sequence logos.

Note: Sequence alignment is a fundamental concept in bioinformatics, and there is a tremendous amount of research and development on the topic. It's hard to imagine doing any amount of serious DNA sequence analysis responsibly without a basic understanding of sequence alignment. I would recommend searching Google and/or Google Scholar for intro-level reviews of the topic.


Error with BWA Mem input having multiple fastq files using cat and process substitution Running BWA Mem on a number of paired end fastq files using process substitution on the inputs results in this error:

[mem_sam_pe] paired reads have different names: "HWI-D00360:5:H814YADXX:1:2111:20284:33068", "HWI-D00360:5:H814YADXX:1:1213:17580:34847"

The command:

bwa mem -M ref.fa "<(cat U0a_CGATGT_L001_R1_004.fastq U0a_CGATGT_L001_R1_005.fastq U0a_CGATGT_L002_R1_001.fastq U0a_CGATGT_L002_R1_002.fastq U0a_CGATGT_L002_R1_003.fastq U0a_CGATGT_L001_R1_001.fastq U0a_CGATGT_L001_R1_002.fastq U0a_CGATGT_L002_R1_004.fastq U0a_CGATGT_L001_R1_003.fastq U0a_CGATGT_L002_R1_005.fastq)" "<(cat U0a_CGATGT_L001_R2_003.fastq U0a_CGATGT_L002_R2_005.fastq U0a_CGATGT_L001_R2_004.fastq U0a_CGATGT_L001_R2_005.fastq U0a_CGATGT_L002_R2_001.fastq U0a_CGATGT_L002_R2_002.fastq U0a_CGATGT_L001_R2_001.fastq U0a_CGATGT_L002_R2_003.fastq U0a_CGATGT_L001_R2_002.fastq U0a_CGATGT_L002_R2_004.fastq)" > out.sam

However the read (there were several) that threw this error does seem to be properly paired:

grep 'HWI-D00360:5:H814YADXX:1:2111:20284:33068' *.fastq
U0a_CGATGT_L001_R1_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 1:N:0:CGATGT
U0a_CGATGT_L001_R2_004.fastq:@HWI-D00360:5:H814YADXX:1:2111:20284:33068 2:N:0:CGATGT

I can see that the concatenation of the input files is happening in a different order, so this is the next thing I'm going to try to fix, but has anyone else run into this?

(I probably just answered my own question but I'll post it anyway, that may or may not be the issue)

You are correct that the concatenation of input read files will confuse BWA. BWA depends on the input reads being in the same order in R1 and R2 files for paired-end alignment.

From bwa docs: If mates.fq is present, this command assumes the i-th read in reads.fq and the i-th read in mates.fq constitute a read pair. BWA will not search through FASTQ files trying to find matching mate IDs.

You will notice that your initial error message is reporting two different read IDs that are improperly paired. Presumably the mispaired read corresponds to the first read in the out-of-order file.


Weird phylogenetic tree I am required to create a phylogenetic tree(s) of my genes of interest and am currently blasting (tblastn) my genes of interest from S. mediterranea against model systems before expanding my search.
Problem I get something like the attached which does not look like right. What am I doing wrong? I'm aligning via MUSCLE, and running Maximum Likliehood Protein models and consequentially, using the correct model (JTT, for instance).
Greatly appreciate any help in advance.

Anyway the chances are this is gene duplication. It will be a housekeeping type gene or else histone that sort of thing. You simply map the different genes onto the various metazoan classes. In your case nematodes and arthropods are basal metazoan, so their is correct given the clade they form part of all all representatives is that gene.

There is nothing weird you just need to read the deep phylogeneticists work like Tom Caveliar-Smith and Sandy Balfour, possibly Blair-Hedges

There are alot of eukaryotes prior to the invertebrates, just look at the Cambrian explosion


How to plot genomic.fna fasta file in R using Gviz? I downloaded the genome of Staph aureus as DNA sequence from NCBI. I would like to visualize it using Gviz. Is Gviz the right tool to visualize genomes? What are the necessary steps from the sequence file to an object that can be used in Gviz? I searched the user guide and the reference manual. It is stated that it is possible to use fasta file, but not how to use them.

We will generate lots of genomes that we will have to annotate and compare. S aureus is not S aureus. There is a spectrum of phenotypes, they all have a different genomes. And we will also not be limited to SA. The question here is not how to annotate the genome, but rather the steps necessary to visualize something and get going. It is my first steps with genomic data, I don't even know what makes sense.

Here is the head of the file:

>NC_007795.1 Staphylococcus aureus subsp. aureus NCTC 8325 chromosome, complete genome
CGATTAAAGATAGAAATACACGATGCGAGCAATCAAATTTCATAACATCACCATGAGTTTGGTCCGAAGCATGAGTGTTT
ACAATGTTCGAACACCTTATACAGTTCTTATACATACTTTATAAATTATTTCCCAAACTGTTTTGATACACTCACTAACA
GATACTCTATAGAAGGAAAAGTTATCCACTTATGCACATTTATAGTTTTCAGAATTGTGGATAATTAGAAATTACACACA
AAGTTATACTATTTTTAGCAACATATTCACAGGTATTTGACATATAGAGAACTGAAAAAGTATAATTGTGTGGATAAGTC
GTCCAACTCATGATTTTATAAGGATTTATTTATTGATTTTTACATAAAAATACTGTGCATAACTAATAAGCAAGATAAAG
TTATCCACCGATTGTTATTAACTTGTGGATAATTATTAACATGGTGTGTTTAGAAGTTATCCACGGCTGTTATTTTTGTG
TATAACTTAAAAATTTAAGAAAGATGGAGTAAATTTATGTCGGAAAAAGAAATTTGGGAAAAAGTGCTTGAAATTGCTCA
AGAAAAATTATCAGCTGTAAGTTACTCAACTTTCCTAAAAGATACTGAGCTTTACACGATTAAAGATGGTGAAGCTATCG
TATTATCGAGTATTCCTTTTAATGCAAATTGGTTAAATCAACAATATGCTGAAATTATCCAAGCAATCTTATTTGATGTT

I am not familiar with Gviz, but I am familiar with a variety of tools for visualizing genome features. At the scale of entire microbial genomes, the primary sequence itself is an overwhelming amount of information and can't really be visualized in a meaningful way by itself. I mean, you could color the nucleotides, but few patterns in DNA can be discerned by an unassisted human eye.

What we're typically interested in are features encoded in the genome sequence: protein-coding genes, non-coding RNAs, transposable elements, regulatory elements, binding sites, epigenetic modifications, and so on. A particular scientist is usually only interested in a subset of these features, based on what particular questions they are investigating. In your comment you say you "would like to annotate protein coding regions and compare them to other genomes", so it seems like we're on the same page here.

Visualization tools often display the genome linearly (even, for example, if it's a circular bacterial chromosome), zoomed out so that the nucleotide sequence isn't visible, and overlayed with boxes, arrows, or other glyphs representing the genes or other features of interest. Here are some examples:

This question thread contains images generated with the GenomeTools AnnotationSketch program (gt sketch command).
The UCSC genome browser which provides web-based access to dynamic displays of the human genome and many other mammalian genomes.
The IGV browser can be downloaded and easily run on your personal laptop or desktop computer.
I'm not familiar with Gviz, but it definitely looks like it belongs in this class of tools.

In each of these cases, the visualization software is not doing the work to determine the positions of the genes and other genome features. This was done independently using tools designed for this task, storing the information in an annotation file in GFF3, GTF, or BED format. The visualization tools need both the genome sequence (in Fasta format) and one or more annotation files (in GFF3, GTF, or BED format) to make meaningful visualizations.


Imputing missing genotypes from separate genotyping panels What is the current standard for imputing missing genotypes between two genotyping panels?  I have two populations genotyped using two different panels (A & B), and I would like to impute all the genotypes in population B for those positions on used in panel A.

I've read the examples for impute2, and I think the closest thing to what I am looking for is this example, "Imputation with one unphased reference panel".

Simply put, I want to provide a list of SNPs, some variant file for population B, and haplotype information from 1,000 Genomes and get imputed genotypes for each SNP in the list.  Is impute2 the state of the art for this?

Given that you mention wanting to use 1000 Genomes as a reference panel for imputing genotypes into your two SNP chip panels, I am going to assume that you are working with human data.

In that case there are several options you can go with:

If your two panels are of European descent, then you are probably best off using the HRC reference panel together with a fast genotype imputation tool such as Beagle 4.1 to impute genotypes in each of your two SNP chip panels separately.
If your panels are not of European descent, then you will likely want to use the 1000 Genomes phase 3 reference panel with Beagle 4.1, Impute2, or Minimac3.

In either case, there are two phasing services available that will do much of the heavy lifting for you 1,2.

The second Wellcome Trust Case-Control Consortium paper performed a cross-imputation analysis as you describe.  I don't see many studies using multiple SNP chip panels.  You will need to take care in your analysis that you are not hit by batch effects from using two different SNP chip panels.

Also, none of these methods will work if the region you are imputing into has too few variants.  I'm not sure what the minimum number of variants is, but if you are using a whole genome genotyping panel of at least 500k SNPs, then you should be ok if you impute a whole chromosome at a time.


Detecting structural variants with MinION data Working on various cancers I have an interest in detecting structural variation (SV) in human, we've successfully used various tools like Pindel, SVDetect, Manta, and LUMPY, to name a few for detecting SVs in illumina short-read sequencing.  I curious if anyone has successfully used ONTs MinION sequencer for detecting SV, as there are many cases where longer reads would be beneficial for SV detection especially where events are long and/or occur in repetitive regions.  Has anyone tried or had success with tools previously designed for PacBio data such as Sniffles?

Supplementary question, obviously the MinION has not quite got the throughput of say the PacBio SEQUEL so can't generate low-pass coverage for the whole human genome as easily, has anyone got experience trying to generate reads targeting specific areas of the genome in areas that are known for translations, long inversions, duplications, etc?

There was a Structural Variant breakout session at the London Calling conference this year. Unfortunately I didn't attend that session, but MinION community members have access to Constance Donnell's summary of that:

https://community.nanoporetech.com/posts/breakout-structural-varia

Here are my attempts at grabbing non-creative chunks from those notes:

Professor Wigard Kloosterman has been developing a bioinformatics pipeline called NanoSV for mapping genomic structural variants in patients with congenital abnormalities
Tomas Sesani at the University of Utah studied the precise copy number of duplicated genes in individual virus genes over time, and tracked the SNP within copy number variable regions of the genome [no public tool mentioned]
Dr. Sudha Rao from Genotypic Technology in India used the MinION sequencer with Sniffles for structural variation calling

ONT have said that the sessions/talks from London Calling 2017 will all be made publicly available at some point in the future.


Separate boxplots for multiple violin plot I am using the following function from seurat package to generate multiple violon plots and I am interested in adding box plots to them but it doesn't work when I have plotted different data at once.
Is there a way to solve it ?

For example, this works:

library(Seurat)
VlnPlot(object = pbmc_small, features.plot = 'PC1') +
  geom_boxplot()

But this will simply lead into an empty box on top of my plots:

VlnPlot(object = pbmc_small, features.plot = c('PC1', 'PC2')) +
  geom_boxplot()

The vioplot package comes built in with boxplots. You can download it from CRAN or there are more features (including formula input and separate colours) in the development version on GitHub: https://github.com/TomKellyGenetics/vioplot

devtools::install_github(“TomKellyGenetics/vioplot”)
library(“vioplot”)
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident)

To plot several PCs:

library(“RColorBrewer”)
colours <- brewer.pal(length(unique(pbmc_small@meta.data$ident)), “Set3”)
par(mfrow = c(1, 2))
vioplot(pbmc_small@dr@pca[,1]~pbmc_small@meta.data$ident, main = “PC1”, xlab = “Identity”, col = colours)
vioplot(pbmc_small@dr@pca[,2]~pbmc_small@meta.data$ident , main = “PC2”, xlab = “Identity”, col = colours)
par(mfrow = c(1, 1))


How can I convert SCF trace files to ABI files? I have a large collection of sequence trace files which are in SCF file format (binary files that in addition to the string of DNA bases also contain the electropherogram of the sample and quality information about the base calls).

I have an existing BioPython pipeline that accepts files in AB1 format (pdf), so I need to convert my SCF files to AB1 files.

I found a question on Biostars that asks for the opposite: converting from SCF to AB1.

The accepted answer there was to use convert_trace from the Staden package.

./convert_trace -out_format scf < trace.ab1 > trace.scf

So I tried the opposite:

./convert_trace -out_format abi < trace.scf > trace.ab1

This did indeed produce a file with *.ab1 extension, but upon importing it to my BioPython pipeline, I got the error

File should start ABIF, not '\xaeZTR'

So apparently, convert_trace converted my files to ZTR format, even altough I specified -out_format abi.

Closer reading of the manual states that [convert_exe] can write CTF, EXP, PLN, SCF and ZTR formats. (so indeed no ABI), but it confuses me why ABI is still a valid output format.

So, my question is, how do I convert SCF files into ABI files?

Apparently convert_trace does not do a good parameter checking and silently sets ZTR as default if it does not recognise a valid output format.

I don't know of any freely available command-line converters in that specific direction. Even at the height of Sanger sequencing (late 90s, early 2000s), pipelines usually tried to reduce complexity and footprint by converting AB1 to anything else, often SCF (ZTR came too late and never got enough traction).

That being said, I dug out an SCF I had from 1997, loaded it into Geneious, and re-exported it as .ab1. I was a bit astonished to see that this worked. So Geneious (or similar packages from other companies) may be your best bet for a quick, one-off thing.

Depending on your needs, you might want to think about other ways: either expand BioPython to read SCFs, or maybe to expand the Staden io_lib to be able to write .ab1 (James still maintains that package and is extremely helpful).

Good luck.


PLINK - Transposed BED file? I'm currently working with the PLINK file system to store large amounts of genotype data.

The plain format consists of three files, two files for storing phenotype and marker information as well as one file (.PED) containing the genotype matrix.

There exists a binary version of the .PED file (.BED), which stores the data in a number of consecutive byte blocks, with each block encoding the genotype of all samples for a specific marker.

While there exists a transposed version of the .PED file (.TPED) which is also supported by PLINK, I couldn't find out whether there also exists a transposed version of the binary .BED file.

Does it exist at all? Are there other file formats to compactly store a genotype matrix for fast sample / marker extraction?

Yes, this exists, and can be efficiently generated by plink 2.0's --export ind-major-bed command.  (The third byte is 0 instead of 1, and the specification is otherwise identical to that of regular plink .bed files with samples and variants swapped.)


Sequence alignment using Markov Model I am learning about applying Markov model to sequence alignment. The prof says that the transition probabilities from a gap-residue alignment to a residue-gap alignment and vice versa are both 0. Is there any biological/mathematical reason behind this statement? Why are the (X,Y) and (Y,X) cell 0? This is a lecture slide of lecture 1, week 4 of the "Bioinformatics: Introdcution and Methods" course on coursera.

If I understand your question correctly, then I think for case of pairwise alignment, there is a simple explanation.

I believe the key insight is that: a mismatch should always score better than a gap.*

This follows biologically since the insertion/deletion (indel) rate is roughly 1/10th that of the substitution rate (i.e. the occurrence of single nucleotide changes), at least in vertebrates. (This varies across the tree of life but I think the substitution rate virtually always exceeds the indel rate.)

To understand why this matters, consider an example:

ATG-AG
ATGT-G

This is an 'impossible alignment' under the probabilities you gave since here we have a transition from a gap-residue alignment to a residue-gap.

However, under our assumption that mismatches are more likely biologically than indels, the correct alignment should be:

ATGAG
ATGTG

Indeed, the latter does look like a better alignment.

This also follows for more complex examples, so this:

ATG--AAG
ATGTT-AG

Becomes this:

ATG-AAG
ATGTTAG

(Or this:

ATGA-AG
ATGTTAG

)

* Strictly, I mean a substitution should score better than an indel (with the associated gap opening and extension penalties). In fact, for the assumption to always be true, a run of mismatches should still score worse than a single indel. This may not always be a correct assumption, consider this example below, is the true alignment case 1) or 2) or something else? Or is in fact a global alignment bad here and this should be split into 2 local alignments? Is there a likely biological mutational event that could explain this? I ask these questions just to point out it is not black-and-white, I don't have clear answers

1)

CGTACGTAGAGGAATGCCCCCCCCC--------AGCAACGTAGCAT
CGTACGTAGAGGAATG---------TTTTTTTTAGCAACGTAGCAT

2)

CGTACGTAGAGGAATGCCCCCCCCCAGCAACGTAGCAT
CGTACGTAGAGGAATGTTTTTTTT-AGCAACGTAGCAT


What is the difference between samtools, bamtools, picard, sambamba and biobambam? After some google searches, I found multiple tools with overlapping functionality for viewing, merging, pileuping, etc. I have not got time to try these tools, so will just see if anyone already know the answer: what is the difference between them? Performance? Features? Or something else? Which one is generally preferred? Samtools?

The obvious answer is that different people wrote them. It's fairly common in bioinformatics for people with a computer science background to get frustrated with existing tools and create their own alternative tool (rather than improving an existing tool). Over time, tools with similar initial aims will have popular functionality implemented in them (and eventually have bugs fixed), such that it matters less which particular tool is used for common methods.

Here's my impression of the tools:

samtools -- originally written by Heng Li (who also wrote BWA). The people who now work on samtools also maintain the alignment file format specification for SAM, BAM, and CRAM, so any new file format features are likely to be implemented in samtools first.
bamtools -- this looks like it was written by Derek Barnett, Erik Garrison, Gabor Marth, Michael Stromberg to mirror the samtools toolkit, but using C++ instead of C
picard -- Java tools written by the Broad Institute for manipulating BAM/SAM files. Being written in Java makes it easier to port to other operating systems, so it may work better on Windows systems. I'm more familiar with picard being used at a filtering level (e.g. removing PCR duplicates), and for statistical analysis, but it links in with the Java HTS library from samtools, so probably shares a lot of the functionality.
sambamba -- a GPL2-licensed toolkit written in the D programming language (presumably by Artem Tarasov and Pjotr Prins). I haven't used it (and don't know people who have used it), but the github page suggests "For almost 5 years the main advantage over samtools was parallelized BAM reading. Finally in March 2017 samtools 1.4 was released, reaching parity on this."
biobambam -- written by German Tischler in C++. I also have no experience with this toolkit. This seems to have some multithreading capability, but is otherwise similar to other toolkits.


Which tools for differential expression analysis in scRNA-Seq? I am starting to run analysis for differential expression in scRNA-Seq. Which tools are available for this kind of analysis?

Can tools for bulk RNA-Seq like DESeq be used for scRNA-Seq?

According to Jaakkola et al. (ref. 1):

  Based on our comparisons, DESeq and Limma without any modifications are not suitable for scRNA-seq data analysis, and yet, they have performed well in the context of bulk RNA-seq data.

They also provided a summary of the comparison between 5 tools (2 designed for Single-Cell, 2 designed for bulk, 1 general):

In my case I ended up using MAST (ref. 2, github).

Jaakkola, M. K., Seyednasrollah, F., Mehmood, A. & Elo, L. L. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics 18, 735–743 (2017).
Finak, G. et al. MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 1183 (2015).


How to share a collection of sbol designs in SynBioHub We currently work on a project that will generate a collection of recombinant plasmids targeting over 3500 genes of the cyanobacterium Synechocystis sp. PCC 6803. We would like to share these designs in SynBioHub but we are not sure how best to go about it. The plasmids have been designed and now they are starting to be produced in batches, checked and finally transformed into Synechocystis to induce a single gene deletion by homologous recombination.
During the project, we will have:

Stage 1 – the initial plasmids designs, waiting for production. They
are worth sharing as the community could influence the production
queue.
Stage 2 – the assembled plasmids. In some cases they will have
different structure than initial designs from stage 1. If nothing
else their description should reflect the status change to
'produced'.
Stage 3 – plasmids are verified by sequencing (also to
check the barcodes are correct) and sequences attached.
Stage 4 – plasmids will be transformed into Synechocystis and checked by PCR to
see if i) the target gene has been knocked out and ii) if the strain
is segregated (i.e. all copies of the genome in the strain contain
the gene deletion).  This information should also be recorded with
the description.

Our initial plan was:

Create a collection with initial designs (component definitions) Stage 1
As stage 2 progresses add new versions of udpated designs
Stage 3, 4 update plasmids descriptions (synbio desc or notes) and attach additional files (e.g. from sequencing).

However, our current understanding is the public collections are frozen so stage 2,3,4 cannot be accommodated.
We are seeking advice how to plan for sharing this resource and some clarifications.

Is there a way of making collection readable by others without freezing its content?
What operations are possible on public component definitions, can the free text descriptions/notes be updated or file attached?
Can new component definitions or their versions be added to a published collection
For 2) 3) who can perform these operations, the collection’s creator or SynBio curators?

Looking forward to input from SynBioHub community

There is a balancing act between allowing all edits, which defeats the purpose of having URIs point to static objects, and allowing no edits, which defeats the purpose of providing some updated information about an object.  SynBioHub resolves this issue by allowing a limited amount of edits after publication.  To answer your questions specifically:

Yes, you can share objects by generating a share link OR adding another owner to an object.  Sharing can be done by the original owner of the object.
There are some mutable fields that any owner of the object can change, these include the description, notes, source, and references in the Details section.  You can also add attachments.
Yes, when you publish (Make Public) a new object or collection, you have the option of making it public into an existing collection OR a new collection.  A person with curator privileges on the SynBioHub instance you are working on can perform this operation.
Answered above.


Doing plot with this data I have normalized expression values of some genes in some patients and the how many days they have survived from their diagnosis with cancer like below

  EDITED

I have taken mean of the expression values for each patient and I have how many days each patient survived since diagnosis as below

>

I done so but nothing I got

> fit <- survfit(Surv(time) ~ gene, data = km)
> print(fit)
Call: survfit(formula = Surv(time) ~ gene, data = km)

                  n events median 0.95LCL 0.95UCL
gene=-0.333140816 1      1    859      NA      NA
gene=-0.307846735 1      1    347      NA      NA
gene=-0.303559694 1      1   1339      NA      NA
gene=-0.290518776 1      1     61      NA      NA

The question is if these genes increase survival time or not but likely I am doing wrong

  By your help

> fit <- coxph(Surv(time) ~ gene, data = km)
> print(fit)
Call:
coxph(formula = Surv(time) ~ gene, data = km)

       coef exp(coef) se(coef)     z     p
gene 0.8664    2.3783   1.4520 0.597 0.551

Likelihood ratio test=0.35  on 1 df, p=0.5543
n= 69, number of events= 69

For plotting I guess I do need another information because

 library("survminer")

> ggsurvplot(fit, data = km)
Error in ggsurvplot(fit, data = a) : object 'ggsurv' not found

  EDITED

I  divided patients based on the median of survival days to Up and Down category like below but still I am failing to visualize that

 fit <- survfit(Surv(Time, Status) ~ Gene, data = km)

> ggsurvplot(fit, data = km)
Error in data.frame(..., check.names = FALSE) :
  arguments imply differing number of rows: 70, 0, 140

Any help please?

Thanks

In general, survival analysis can be said to be composed of two steps; Cox regression, with which you calculate the "hazard ratio" based on your variables, and a "Kaplan-Meier (KM) estimate", which is used to visuazlize the data.

Here is a nice tutorial for doing survival analysis with the survival and survminer packages. The latter includes ggplot2 kind visualization of the KM curves.

For the Cox regression, you can use continuous(*) or discrete variables, in your case these would be the normalized expression value of your gene and categories of your gene respectively. You can categorize the expression level of a gene as low vs high, or low vs medium vs high, based on different quantiles.  Another approach to categorize your expression data would be to find an "optimal" cut point with maximally selected rank statistics; from the link:

  This is an outcome-oriented methods providing a value of a cutpoint
  that correspond to the most significant relation with outcome (here,
  survival)

Basically you search for the cut-point where your data is split into two such that the two resulting parts of the data show the most separation/difference in terms of survival.

* Can a Cox Proportional-Hazards Model be built only with continuous predictors?


GATK 4.2 using IndexFeatureFile with multiple files I am trying to use GATK IndexFeatureFile to create ".vcf.gz.tbi" file. But I have 480 files.
For the pilot study I used:
!./gatk IndexFeatureFile -I HG00096.chrom22.bam.10061.vcf.gz -O HG00096.chrom22.bam.10061.vcf.gz.tbi

Doing this procedure for 480 files is not seem very practical. Is there any way to use all file as an input and it takes all files and gives ".tbi" files as an output?
I thought
!./gatk IndexFeatureFile -I *.vcf.gz -O *.vcf.gz.tbi

or maybe not using output argument (just input). I did not worked for me.

You could run loop over the files in your shell, running a separate command for each file. For example, in Bash (or another Bourne shell):
for file in *.vcf.gz; do ./gatk IndexFeatureFile -I ${file} -O ${file}.tbi; done


Single-cell sequencing dataset has too many barcodes I am analyzing a single-cell sequencing dataset from the website 10xgenomics, with 2000 cells. It is a BAM file and I am trying to obtain the individual cells per sample. I used the command

samtools view mtdnaAsorted.bam | grep CB:Z: | sed 's/.CB:Z:([ACGT]).*/\1/' | sort | uniq -c > reads_per_barcode_mtdnaA

to get a list of the individual bar codes and their reads, but there are 100,000 unique barcodes even with only 2000 cells. Does anyone know what could be the problem?

I know that there could be noise, but this is a problem for me because I'm trying to analyze coverage per cell, and if there are a lot of cells that contribute towards the total number of reads, the net coverage per cell is going to be really low

The BAM file you downloaded has all of the detected 10X barcodes instead of those that represent cells.  In a given 10X experiment the number of input barcodes vastly outnumbers the input cell count.  After the reaction some of barcoded beads are invariably sequenced because they are encapsulated with ambient RNA, dead cells etc. After processing the BAM file Cell Ranger infers which cells likely represent cells and generates gene expression matrices from these files.

You would need to download one of the output files that contains the list of passed barcodes such as the Gene / cell matrix (filtered) and then use the list of passed barocode sequences to filter your BAM file. Alternatively you could approximately select the same cells by selecting the top 2,000 barcodes with the highest numbers of reads.


different gene expression in the monocle As it is shown in the picture that we can get the pseudotime plot by using Monocle. I was wondering whether we can use some function like the Seurat FindMarker to find out the differential expression genes between the Double-KO and WT in the  bottom right corner?

Seurat FindMarker Doc
Yes based upon the R documentation here it seems possible.

This  Differential expression article  seems to show this is how.

library(Seurat)
pbmc <- readRDS(file = "../data/pbmc3k_final.rds")
# list options for groups to perform differential expression on
levels(pbmc)
# Find differentially expressed features between CD14+ and FCGR3A+ Monocytes
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono")
# view results
head(monocyte.de.markers)
# Find differentially expressed features between CD14+ Monocytes and all other cells, only
# search for positive markers
monocyte.de.markers <- FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = NULL, only.pos = TRUE)
# view results
head(monocyte.de.markers)
# Pre-filter features that are detected at <50% frequency in either CD14+ Monocytes or FCGR3A+
# Monocytes
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", min.pct = 0.5))
#...
# Test for DE features using the MAST package
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "MAST"))
head(FindMarkers(pbmc, ident.1 = "CD14+ Mono", ident.2 = "FCGR3A+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Maybe something like this would work for you.  Monocle export

  #replace the monocle_cds with your monocle
  seurat <-exportCDS(monocle_cds, export_to = c("Seurat", "Scater"))
  #This bellow will list the options for ident.1 and ident.2
  levels(seurat)
  # insert name from levels(seurat) command in parentheses
  head(FindMarkers(seurat, ident.1 = "Double-KO+ Mono", ident.2 = "WT+ Mono", test.use = "DESeq2", max.cells.per.ident = 50))

Hopefully you are able to separate the two groups Double-KO and WT into ident.1 and ident.2 so you can use the function on them. Good luck. As long as the exportCDS function works on your data type you input for cds.


How to quantile normalization on RNA seq counts I have a read count data (RNAseq) and want to perform quantile normalization. Could you please help me how to do it. I tried some scripts in R but it didn't work. I want the result output in matrix form.

  gene_id   SRR896664   SRR896663   SRR896665
  ENSG00000000003   46106   36353   40614
  ENSG00000000005   198 399 1200
  ENSG00000000419   40364   37769   40849
  ENSG00000000457   18924   16211   16057
  ENSG00000000460   31040   28888   29901
  ENSG00000000938   200 0   394
  ENSG00000000971   14935   14353   12522

The script which i tried is

  data <- read.csv("testquantile.csv",header=T)
  head(data)
  rownames(data) <- data[,1]
  data_mat <- data.matrix(data[,-1])
  head(data_mat)
  data_norm <- normalize.quantiles(m, copy = TRUE)


running a samtools command for multiple bam files from 1000 genomes project I want to compute the depth of coverage only for specific intervals in phase 3, 1000 genomes project.I have not worked with 1000 genomes project before, so a bit unfamiliar with it. I do not want to download all of the bam files for the entire genomic region, just those intervals which I am interested in them.
I guess I can use this (I found from biostars forum)

samtools -bu view 'http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam' "2:1000000-2000000" | (...)

However, I am just a bit confused. First of all, when I try the code above directly from the command line, I get an error message "No such file or directory"! Am I doing something wrong? or should I provide a specific path?
Second, I want to run  it for the mapped files for all individuals: *.bam”. through a loop.

Can anyone help me to solve my problem?

The command you're looking for is:

samtools view -c http://ftp.1000genomes.ebi.ac.uk/vol1/ftp/phase3/data/HG01375/alignment/HG01375.mapped.ILLUMINA.bwa.CLM.low_coverage.20120522.bam 2:1000000-2000000

For a series of intervals, you'd be better off scripting something with pysam:

import pysam
bam = pysam.AlignmentFile("some file.bam")
regions = open("your regions.bed")
for region in regions:
    chrom, start, end = region.split("\t")
     bam.count(chrom, int(start), int(end))  # Do something with the result

You might see if there are bigWig files available, they'd be faster to process.


Snakemake wildcard constraints reject too much I have the following two rules:
rule extract:
    input: "data/{dir}/{file}.gz"
    output: "data/{dir}/{file,^.*(?<!\.gz)$}"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k {wildcards.file}.gz"

rule extract_dot:
    input: "data/{dir}/wtdbg2.3.dot.gz"
    output: "data/{dir}/wtdbg2.3.dot"
    conda: "config/conda-extract-env.yml"
    shell: "cd 'data/{wildcards.dir}'; gunzip -k wtdbg2.3.dot.gz"

The second one is actually a hack because the first does not work on files named wtdbg2.3.dot.
I have also tried giving the constraint using the wildcard_constraints member of the rule, and prefixed the string with an r or not.
Whatever I do, as soon as I uncomment the extract_dot rule, snakemake complains about the missing input file data/C.elegans/wtdbg2.3.dot in another rule. Note that there is a rule with output data/C.elegans/wtdbg2.3.dot.gz, and when changing this rule to output data/C.elegans/wtdbg2.3.dot instead, everything works as well.
So I conclude that the problem is in the regex. I have tested the regex on regex101.com, and it accepts wtdbg2.3.dot. But snakemake seems not to. Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

Does anyone know if there is anything special about the snakemake regexes? Should they not be normal python regexes?

No, nothing special.  The Snakemake documentation just links to the Python re module documentation.
I think this is something strange in how the anchors ^ and $ behave when matching wildcards, and not anything else in your particular regular expression.  Check out the pattern-matching with this dummy Snakefile:
rule all:
    input: ["data/subdir/example", "data/example", "example"]

rule anchored_example_of_subdir:
    """Can't add either ^ or $ or this will fail"""
    output: "data/{subdir,subdir}/example"

rule anchored_example_in_subdir:
    """Can't add ^ or this wil fail"""
    output: "data/{file,example$}"

rule anchored_example:
    """This works with both ^ and $"""
output: "{file,^example$}"

(There are no inputs or actions so it's just a matter of what filenames match what rules, and the regular expressions do nothing but match exact strings.)  This "works" as is (snakemake -nrp shows what you expect) but if you add in anchors as described in the docstrings above the upper two rules no longer match the expected filenames.  It's as if (and maybe is?) under the hood Snakemake is trying to apply the expression to the entire path.
Can you just drop your ^ and $?  The snakemake documentation implies that it's using the equivalent of re.match rather than re.search already, with wording like "Here we could restrict the wildcard dataset to consist of digits only using \d+ as the corresponding regular expression" so I don't think it's necessary anyway.

Edit: I think I found the issue in the snakemake source.  It munges your constraints together into one larger regular expression for the entire file path, so ^ and $ get embedded inside the resulting expression.  This looks like something that should at least be clarified in the documentation.