How to make a robust amino acid phylogeny with adequate robustness support This is the second part of this thread:  Duplicated genes and genome fragmentation

As @Michael G. suggested, I did a pseudo-annotation of sequences 5-2 and 5-3 (i.e., I joined both sequences into one) and I did the phylogenetic analysis again. Here is the tree that I obtained:

Both of my sequences are shown in a red rectangle. The length of the branch seems to have decreased a bit. Also, I looked at the branches near sequence 5-1. One of them is the original query sequence from species A and the other one is a hypothetical protein from the same species (probably related to the query sequence, but I couldn't find any domains because it is only a few amino acids long). I still need to find the root of the phylogeny but I am not sure how to do this...

To finish you need bootstrap support using a model-based tree building algorhithm, I use RAxML, specifically raxmlHPC (easily downloadable and compiles on Linux and OSX). The codes are quite complicated and I've given them below.

A robust maximum likelihood tree is,

raxmlHPC -m PROTGAMMAILG -n output.tre -o outgroup-p 10000 -s alignment_file.phy

To make a bootstrap tree,

raxmlHPC -b 459 -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -s /path_to_alignment.phy -N 500

Make a consensus tree of the bootstraps,

raxmlHPC -f b -m PROTGAMMAILG -n output.tre -o outgroup -p 10000 -t /path_to_master_tree(best_tree) -z path_to_bootfile

You require bootstrap support >80% and please repeat this with and without 5-2/5-3 (it still looks long)

You need a fairly robust desktop, dedicated for around 24 hours to complete one calculation. You have enough results to draw a conclusion.

Viewing the tree, FigTree is easy.

Root

Check the work of Cavalier-Smith, he heavily focused on basal protists to find a species which evolved, using external evidence, before any other member of your group. One of the ingroups may have sufficient evidence to say it evolved before the rest and that is your root, which you select as an outgroup.


Mapping protein refseq to Gene ID I have a protein refseq (eg, NP_000029). How can I get the corresponding gene ID (ag, APC) from NCBI using an R package?

You can do this with the reutils package, which provides an API to NCBI's E-utilities. Here's an example for your specific question:
install.packages("reutils")
library(reutils)

# Get universal identifier
uid <- esearch("NP_000029", db = "gene")
# Fetch summary
sm <- esummary(uid, db = "gene")
# Extract specific tag
sm$xmlValue("//Name")


How BLOSUM Matrix is constructed and calculated I would like to ask how BLOSUM matrix is constructed and calculated ? I read the wikipedia and I still do not understand it. I do not understand the mathematical calculations as I have low knowledge about logarithms.

Maybe could someone give a simple example and then a real-scenario example of how the BLOSUM matrix can be used and calculated, maybe in relation to the Blast algorithm where it uses the BLOSUM scoring matrix to determine high scoring words for each word in the query sequence ?

Also, why are they different types of BLOSUM matrices like BLOSUM62 and so on ?

Sean Eddy wrote an excellent explanation of the BLOSUM Matrix in his Nature Biotechnology paper Where did the BLOSUM62 alignment score matrix come from? He also gives real examples of how one would calculate the score for different amino acid pairs and where the numbers come from.  Further references given at the end of the paper might also be useful reading, specifically: Amino acid substitution matrices from an information theoretic perspective, Methods for assessing the statistical significance of molecular sequence features by using general scoring schemes, and Amino acid substitution matrices from protein blocks.


Retrieving a list of human genes having GO associations I need a command in R to retrieve all human genes associated with a Gene Ontology entry. I tried to look for it online but did not find it.

Here's an example for the mouse genome:

library(org.Mm.eg.db)
select(org.Mm.eg.db, c("GO:0048406"), c("GENENAME","SYMBOL"), c("GO"))

You get output like:

                                                    GENENAME SYMBOL
1                                     pregnancy zone protein    Pzp
2 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
3 nerve growth factor receptor (TNFR superfamily, member 16)   Ngfr
4                                             neurotrophin 3   Ntf3
5             neurotrophic tyrosine kinase, receptor, type 1  Ntrk1
6            furin (paired basic amino acid cleaving enzyme)  Furin
7              proprotein convertase subtilisin/kexin type 6  Pcsk6
8                                                 sortilin 1  Sort1
9                                      alpha-2-macroglobulin    A2m

Edit: Given your comment that you want instead a list of all genes with a GO term associated with them (for human):

library(org.Hs.eg.db) # Install it from bioconductor yourself
unlist(as.list(org.Hs.egSYMBOL)[mappedkeys(org.Hs.egGO)])

The mappedkey() part gets the index into the bimap of all entries with a valid GO mapping. That's then used to subset the symbol bimap (after converting to a list).


How to perform functional analysis on a gene list in R? From an RNA-seq experiment I have about 17000 gene ids for 2 sample conditions arranged according to their log2 fold changes when compared to a control. I need to annotate these, but I've never done annotation before and am wondering how to do this in R.

I'm interested in up and down regulation of the pathways of the involved genes. I'd love to get enrichment and p-value analysis too. I'm interested in whether or not a package similar to PANTHER exists in R for this analysis. I've been reading how different annotation programs can be better or worse; it seems PANTHER is much better than DAVID, so I was wondering if there is a package that provides PANTHER-like analyses in R.

There seems to be multiple packages available, are there any that stand out as being the best?

I'm primarily interested in human samples and annotated pathways.

N.B., I'm avoiding discussion of "best", since that's more or less impossible to answer.

Your question can actually be divided into two:

What's a good tool for pathways analysis (ideally in R)?
What are good sources of pathway information?

For (1), there are a few different possibilities, but I prefer either roast or camera from limma or goseq (a stand-alone bioconductor package). Of those, I expect what roast or goseq are doing are the most similar to what panther is doing in the sense that they're not competitive. However, since genes do correlate with each other, I prefer competitive tests like that offered by camera.

For (2), your options are basically KEGG or reactome. KEGG is problematic, since you need a license for anything remotely recent (there are some R packages to get around this, but I haven't a clue how kosher they are legally). Given that, reactome might be your best bet.

Having written all of that, the best database out there is what IPA uses. This is a commercial product (I have no affiliation), but you can get a demo license that will work for a couple weeks at least. If you're doing analyses relatively often then it makes sense to spend a bit of cash on a license (maybe share with neighboring groups).

Regardless of what you use, you might also want to use pathview, which can produce easily interpretable plots related to pathway enrichment.


How can I combine a kinship matrix with subset individuals when using rvtests? I've been processing data for a consortium project that uses the rvtests toolkit. Our data analysis process uses a kinship matrix and the calculation of inverse-normalised transformed statistics after covariate adjustment (for 1kG-imputed data from a SNPchip). One of the operations we have been asked to do is to analyse males and females separately for a number of different phenotypes.

Here's an example of the command line that I'm running:

rvtest --inVcf chr${ch}.imputed.poly.vcf.gz \
 --pheno ../scripts/phenotypes.raw.2017-May-01.ped \
 --pheno-name ${traitName} \
 --out rv_results/${traitName}_FEMALE_chr${ch} \
 --kinship kinship_matrix.kinship \
 --meta score,cov[windowSize=500000] \
 --peopleIncludeFile ../scripts/female.iids.txt \
 --covar ../scripts/covar.2017-May-11.ped \
 --covar-name ${covs} --useResidualAsPhenotype --dosage DS

Unfortunately, the kinship matrix doesn't interact well with the --peopleinclude parameter. When I run this, the program complains about unexpected columns in the kinship matrix (starting at one column greater than the number of females in the dataset). I suspect what is happening is that rvtests is assuming that the kinship matrix applies to the sample set after filtering, rather than applying to the total sample set. This sets off a few warning bells in my head, because it means that rvtest is probably ignoring the column/row labels in the kinship matrix that it has generated.

How can I work around this problem?

Adam Locke (a collaborator of this project) suggests that removing covariate information for the unselected individuals (i.e. setting it to NA) works around this problem:

  I believe the problem is that he is using a pre-computed kinship matrix including both males and females, and when using the “—peopleIncludeFile” it can’t properly select the right people from the kinship matrix. I don’t know that others have tried to do analysis this way. An alternative, that I think has worked for others, is to generate phenotype variables that have values only for males or females (so, e.g., BMI, BMI_females, BMI_males). For whatever reason, this appears to correctly only use the and select the samples desired without the kinship matrix.

Presumably a kinship matrix that is generated for just the subset individuals will work as well, although I don't like the idea that rvtests isn't able to do a label lookup on the kinship matrix to find the correct individuals.


data storage requirements for RNA-seq and WGS Supposing upcoming RNA-seq for 3 conditions and 3 replications for each, how much is data storage requirements?

The breakdown would be:

1.5-2 GB per sample for fastq files.
~3GB per BAM file, depending on the aligner

So ~5GB per BAM file, which is about 50GB total. That can vary by maybe 20%, depending on how compressible the results are (I rearrange our fastq files to make them more compressible, so I'm probably underestimating what others would see).

Obviously using Salmon/Kallisto would halve the storage requirements. Using CRAM would also decrease things, though obviously fewer tools support it.


scater: SingleCellExperiment, Error in seq_len(ncol(assay)) I am trying to go through the following scRNA-seq tutorial.

But the line sce <- newSCESet(countData=all.counts) is not working anymore with the most up-to-date version of scater. Now we should use SingleCellExperiment function instead, so when trying to do the following:

sce <- SingleCellExperiment(assays = list(counts = all.counts))

I am getting the error:

Error in seq_len(ncol(assay)) :
  argument must be coercible to non-negative integer
In addition: Warning message:
In seq_len(ncol(assay)) : first element used of 'length.out' argument

Any help would be greatly appreciated.

I think you might want to use the most recent workflow, which should be for the latest version on Bioconductor.


Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype? So if I attempt this question using their method I put the data of genotype and phenotype into a data frame in R then use lm() to do linear regression and summary() to get an F statistic with p value 0.54. I'm not connecting what the F statistic means in this case to the answer.

Question 2: Assuming that we have the following SNP and phenotype data, is the SNP significantly associated with the phenotype? Here, we represent each SNP site as the number of minor alleles on that locus, so 0 and 2 are for major and minor homozygous sites, respectively, and 1 is for the heterozygous sites. We also assume that minor alleles contribute to the phenotype and the effect is additive. In other words, the effect from a minor homozygous site should be twice as large as that from a heterozygous site. You may use any test methods introduced in the chapter. How about permutation tests?

  GENOTYPE: 1, 2, 2, 1, 1, 0, 2, 0, 1, 0

  PHENOTYPE: 0.53, 0.78, 0.81, -0.23, -0.73, 0.81, 0.27, 2.59, 1.84, 0.03

Answer: Take the linear regression with F-test as an example, the p-value without the permutation test is 0.54 and is 0.55 after 1000 permutation (this number could be different with your answer because of the randomness of the permutation test).

Link to where I'm getting this: http://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1002828#pcbi.1002828-Balding1

Your approach is correct - it's not the only one permitted under their conditions, but it would be the standard way given only that information.

F-test compares a "full model" (predicting the phenotype using genotype and an intercept) versus a "reduced model" (predicting the phenotype using only the intercept, i.e. the mean value of phenotype). The p-value tells whether the full model predicts better, i.e. whether adding the genotype variable significantly improves the prediction. In this single-variable case, it should be approximately equal to the p-value for the beta coefficient of the genotype, if that's more familiar.

Anyway, this isn't specific to genetics, so you can check out other online sources - e.g. https://onlinecourses.science.psu.edu/stat501/node/295 for a longer, but accessible explanation of the F-test, and https://stats.stackexchange.com/ for many answers to similar questions.


Is spark widely used in bioinformatics? I learned that GATK 4 is using Spark for parallelization. I googled around, though I am still not quite sure how spark really works and how to use it in practice. Besides GATK 4, are any other bioinformatics tool using spark? Generally, is spark widely used? Is it a necessary skill to learn? Thanks in advance.

The only bioinformatics tool other than GATK4 that I am aware of that uses spark is Hail (a Spark based replacement for Plink). Hail is also supported by researchers at the Broad.

Most places I have worked at have not switched over to Spark.  As such, I don't think it is widely used generally. Therefore, I don't think knowing Spark qualifies as a necessary skill to learn for bioinformatics at this time.

However, I think Spark is a superior paradigm for handling the large amount of data that bioinformaticians routinely deal with, and I think we will soon see the field move towards using it more.  Any bioinformatician would do well to acquaint themselves with Spark and to play around with Hail.


How do you generate read-length vs read-quality plot for long-read sequencing data (e.g., MinION)? How do you generate read-length vs read-quality plot (heat map with histograms in the margin) for long-read sequencing data from the Oxford Nanopore Technologies (ONT) MinION? The MinKNOW software from ONT provides a plot like this during base calling.

This would also be very helpful for PacBio data.

I also wrote a package to create various plots from Oxford Nanopore sequencing data and alignments: NanoPlot. It can be installed through pip (see also the README on Github). In addition to multiple plots also a limited NanoStats output is created (see also NanoStat). Data can be presented using:

A fastq file (optionally compressed)
A bam file
The sequencing_summary.txt file generated by albacore

Using optional flags you can:

Log transform the read lengths
Use aligned reads rather than sequenced reads
Downsample the reads
Set a maximum read length

I've added an example below, plotting log transformed read length vs average read quality (using a kernel density estimate). More examples can be found in the gallery on my blog.

I welcome all feedback and suggestions!


Pipeline for extracting gene from multiple genomes for use in HyPhy selection analyses? I have been trying to obtain some preliminary data from HyPhy selection analyses to inform a larger project. I have obtained a number of assembled mammalian genomes from NCBI with the initial goal of extracting from each of them their corresponding sequence for a specific gene. I wanted to use HyPhy to determine if there are signatures of positive selection in any of the mammalian lineages for that specific gene. However, I have been struggling to make headway with this for some time and am realizing that I may be out of my depth.

I was hoping that people could suggest, or at least point me in the direction of, a pipeline/guide/best practices for how one would go about achieving what I am trying to do. I would like to start from scratch and make sure I am performing my analyses in a proper way, using advice from other users.

As I mentioned, I currently have a collection of assembled genomes from NCBI which I would like to pull specific genes from. My understanding is that I will then need to extract the exons (I imagine there is a way to do this without having to extract the entire gene first), align the sequences from each species while preserving reading frames, and then use that alignment as input in HyPhy's codon models. It is very possible that I am missing out some fundamental steps in the process.

If anybody could offer direction/advice, I would be very appreciative.

Okay well, because you're using HYPHY, you're going to need a gene tree of every gene you intend to analyze. This requires you to understand the homology relationships between the genomes at hand. Depending on your project, you might already have a list of homologs you're trying to analyze, or you can start from scratch.

If starting from scratch, you will need the proteomes of the species you're interested in (that is, the translation of the coding sequences of each gene, or to put it another way, the concatenated exons). Genome releases should typically have a proteome available for download, along with its corresponding coding sequence.

Next, you need some way to get the homolog groups for the genes in your genomes. You can use a pipeline like OrthoFinder, or iterative tree building and trimming, to infer homolog clusters (OrthoFinder calls them orthogroups but they're the same thing - homologous genes for the species of interest including both single copy orthologs and paralogs).

For this step, if you already have a set of genes and you know the copy numbers, or know that they're single copy orthologs, then it is more easy - you can just pull the sequence from each proteome and go from there.

Once you have your homologs, you need to align the amino acids using a tool like MAFFT, and then align the coding sequences to the amino acid alignment to attain a coding sequence alignment. For this you can use a tool like pal2nal, or check out pxaa2cdn from the phyx package written by some of my collaborators.

If the sequences aren't too divergent, and intron-exon structure is conserved, you could possibly also include alignments of the introns to provide more phylogenetic information. If you do use this approach however, you'd need to align the exons and introns separately. Also, the introns would not be included in the selection analysis.

Next, I would typically reinfer the gene tree of each gene using nucleotide models, which depending on the time, level of selective constraint, etc. can give more signal for the phylogenetic relationships. For this you can use a tree inference program like RAxML. With this tree and your coding sequence alignment, you'll be ready to run analyses in HYPHY (aBS-REL, MEME, FUBAR, etc., depending on your requirements).

Note that if you want, you can also do this in a Bayesian framework to integrate over uncertainty in the model parameters, tree topology, etc. Several Bayesian phylogenetic inference tools (e.g. MrBayes, PhyloBayes), allow you to infer putatively positively selected sites at the same time as tree inference.


TruSeq Exome Kit bed file old version I have to run an analysis using Whole Exome Sequencing (WXS) data prepared with the TruSeq Exome Kit (2012).

I need the .bed file for the kit and I found it here.

However, this file was posted on Nov. 16, 2015, while the kit used in my case is from 2012. I searched but I didn't find the bedfile for that specific version online.

Will I run into problems if I use this version, and in this case, where can I retrieve the bedfile for 2012's version?

Usually it is more problematic to work with older versions of data with newer protocols than the other way around.

In this case, I think that the features annotated in the (newer) bed file should be more accurate than the older one. Unless you specifically want to compare your experiment with a previous experiment performed with a bed file from 2012 you'll be fine using the newer one.

However make sure that you are correctly using the right file and that you report (or record) each version of protocol and data used in case some doubts arise.


How to remove degenerate sequence (N4) at the 5' and 3' of the reads? I have to analyse GSE114327. For that, I use GSE instruction based on GSM3139597.

In that data processing section has written:

  FASTQ reads were subjected to adapter removal (cutadapt) per manufacturer's instructions (NEXTflex Small RNA Sequencing Kit v3).

So, I use cutadapt for trimming. Then

  The degenerate sequences (N4) at the 5' and 3' end of the reads were removed

based on the instruction. But I don’t understand that sentence and I don’t know how I can remove degenerate sequence (N4) at the 5' and 3' of the reads.

From the link you posted it looks as though they used UMI tools as per the section which says
'were removed and added to the name of the read as a unique molecular identifier (UMI) using UMI tools'.  The link above has instructions for installation and use.
From this I would guess that this data has Unique Molecular Identifiers (UMIs) for each read.
Looking up the protocol for the nextflex kit used (See Link) confirms this:

  Data Analysis
  The 3' and 5' adapters included in this kit both contain 4 random bases that will appear immediately 5' and 3' to the insert in sequencing data. The presence of these random bases should be considered when choosing an alignment strategy. When using “end-to-end” alignment, we
  recommend processing data in the following manner:
  1. Clip the 3' adapter sequence (TGGAATTCTCGGGTGCCAAGG).
  2. Trim the first and last 4 bases from the adapter-clipped reads.
  3. Perform alignments as normal.

It seems that N4 is a 4 base random barcode at both the 5'& 3' end of the read - these are cut from the read and this metadata is added to the read name so that every read can be uniquely identified.

Edit:

You can install UMI tools via:

conda install -c bioconda -c conda-forge umi_tools

or

pip install umi_tools

Usage instructions can be found at:

https://umi-tools.readthedocs.io/en/latest/reference/extract.html

  Usage:
  For single ended reads, the following reads from stdin and outputs to stdout:

  umi_tools extract --extract-method=string --bc-pattern=[PATTERN] -L extract.log [OPTIONS]

  For paired end reads, the following reads end one from stdin and end two from FASTQIN and outputs end one to stdin and end two to FASTQOUT:

  umi_tools extract --extract-method=string --bc-pattern=[PATTERN] --bc-pattern2=[PATTERN] --read2-in=[FASTQIN] --read2-out=[FASTQOUT] -L extract.log [OPTIONS]


How to trim adapter sequences from GSE65360 in order to map the reads? I am trying to map single-cell chromatin accessibility data from doi:10.1038/nature14590, obtained using scATAC-seq, to the reference genome. Example paired-end reads are here.

The reads contain various adapter sequences, which the authors "trimmed from FASTQs using custom python scripts" [ref]. I am having problems reproducing this trimming. They used Nextera Transposase Adapters

  TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
  GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG

and custom barcode adapters which are listed [here]. An excerpt from this list:

  Custom Barcodes Adapter 1 (index i5):
  AATGATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACTAACCAAGTCGTCGGCAGCGTCAGATGTGTAT
  AATGATACGGCGACCACCGAGATCTACACAGGTTGCCTCGTCGGCAGCGTCAGATGTGTAT
  Custom Barcodes Adapter 2 (index i7):
  CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCCCAACCTGTCTCGTGGGCTCGGAGATGTG
  CAAGCAGAAGACGGCATACGAGATCACCACACGTCTCGTGGGCTCGGAGATGTG

Analyzing the FASTA-files with FastQC showed that the reverse complement of those sequences was over-represented in the data and also confirmed Nextera Transposase Adapters.

I tried to trim the sequences using cutadapt:

cutadapt --minimum-length=20 -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATA -A ATACACATCTGACGCTGCCGACGA -A CTGTCTCTTATA -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz SRR1780164_1.fastq.gz SRR1780164_2.fastq.gz

with

CACATCTCCGAGCCCACGAGAC being the reverse complement of the 2nd fix part of the i7 adapters
ATACACATCTGACGCTGCCGACGA being the reverse complement of the 2nd fix part of the i5 adapters
CTGTCTCTTATA for the Nextera Transposase Adapters

FastQC does not indicate overrepresented sequences in the resulting FASTA-files, however still indicates Nextera Transposase Adapters in the _1 files. On the other hand, after TrimGalore removes CTGTCTCTTATA, Nextra Adapters are not indicated anymore, for some reason.

Furthermore, mapping the reads using bowtie2

bowtie2 -x GRCh38_primary_assembly_index -1 SRR1780164_1_trimmed.fastq.gz -2 SRR1780164_2_trimmed.fastq.gz > SRR1780164.sam

results in a low alignment rate

1352 (2.79%) aligned concordantly exactly 1 time
3420 (7.07%) aligned concordantly >1 times
11.25% overall alignment rate

which gives me the impression I did something wrong during adapter trimming.

Any help would be much appreciated.

--

Edit: As I am not sure about the sequences to use for trimming, I also tried -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA, as found in a suggestion for Nextera Adapters, however it did not remove the sequence over-representation in the -1 file. -a CACATCTCCGAGCCCACGAGAC -a CTGTCTCTTATACACATCTCCGAGCCCACGAGAC -A CTGTCTCTTATACACATCTGACGCTGCCGACGA solves adapter-issues according to FastQC, but also does not give good mapping results.

First the Nextera adapters and the custom barcode adapters overlap each other

Nextera 1                                     TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
Barcode 1   GATACGGCGACCACCGAGATCTACACTAGATCGCTCGTCGGCAGCGTCAGATGTGTAT

Nextera 2                                    GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
Barcode 2    CAAGCAGAAGACGGCATACGAGATACCCAGCAGTCTCGTGGGCTCGGAGATGTG

Looking at some read 1s and read 2s you were correct that the adapters appear at the 3' end of the read reverse complimented (The |> points to the adapter in the sequence).  You can see that CTGTCTCTTATACACA would be the ideal sequence to use for the 3' end of the reads.  Since the Nextera adapter and the barcode overlap you only need to use the one that is adjacent to the insert.

SRR1780164_2.fastq: GCCCGAGGGGACGCTGCCCCACGTGCTCCAGG |> CTGTCTCTTATACACATCTGACGCTGCCGACGATGCACGAAGTG
SRR1780164_1.fastq: GGGCTGGACAGCGGCGGGGACGCGGACCTGCAG |> CTGTCTCTTATACACATCTCCGAGCCCACGAGACACAAAGTGA

Of the 64k reads there are ~11k with the Nextera adapter at the 3' end based on a rough sequence grep

On the 5' end of the read there is some Nextera adapter present but it is not reverse complimented and there are only a handful of reads with the adapter present ~250 total between read 1 and read 2.

SRR1780164_2.fastq: TGGGCTCGGAGATGTGTATAAGAGACAG <| CTCCCCCGCCTCCAGCATCCGGGCGAGGTAGTGCATCGACGCGTCCAC

I used this command:

cutadapt --minimum-length=20 -A CTGTCTCTTATACACA -a CTGTCTCTTATACACA -G GATGTGTATAAGAGACAG
-g GATGTGTATAAGAGACAG -o SRR1780164_1_trimmed.fastq.gz -p SRR1780164_2_trimmed.fastq.gz
SRR1780164_1.fastq SRR1780164_2.fastq

Which gives a trimming result of

   === Summary ===

Total read pairs processed:             64,832
  Read 1 with adapter:                  16,721 (25.8%)
  Read 2 with adapter:                  16,996 (26.2%)
Pairs that were too short:                  46 (0.1%)
Pairs written (passing filters):        64,786 (99.9%)

I think this best represents the library construct and the trimming of the adapter sequences


Designing a lab NGS file database schema I am the resident Bioinfo Geek in a hospital academic lab that routinely employs NGS as well as CyTOF and other large volume data producing technologies. I am sick of our current "protocol" for metadata collection and association with the final products (miriad excel sheets and a couple poorly designed RedCap DBs).

I want to implement a central structured, controlled datastore that will take care of this. I know that the interface to the technicians how will be inputing the data is crucial to its adoption, but this is not the focus of THIS particular question: Does there exist a schema or schema guidelines for this type of database?

I would rather use a model that has been developed by people who know how to do this well. I know of BioSQL but it seems more geared towards full protein/nucleotide records like those found in uniprot or genbank. That is not what we have here. What I want is something similar to the system touched on in this preprint: http://biorxiv.org/content/early/2017/05/10/136358

Alternatively, can anyone provide links to where I might find relevant guidelines or supply personal advice?

The Global Alliance for Genomics and Health has been working on the issue of representing sequencing data and metadata for storage and sharing for quite some time, though with mixed results.  They do offer a model and API for storing NGS data in their GitHub repository, but it can be a bit of a pain to get a high-level view.  I am not sure if any better representation of this exists elsewhere.

I can say from personal experience (having built over a dozen genomic databases), there is no ideal data model and storage best practices.  Genomic data comes in many shapes and sizes, and your needs are going to vary from every other organization, so what works for one bioinformatics group won't necessarily work for you.  The best thing to do is design and implement a model that will cover all of the data types in your workflow and downstream analyses you might do with the data and metadata.


Differences in webtool and standalone Clustal Omega Percent Identity Matrix I'm trying to align sequences using MUSCLE in my own computer. To compute Percent Identity Matrix I'm trying to use Clustal2.1 like it's done in https://www.ebi.ac.uk/Tools/msa/muscle/. Even though I'm doing the same thing, I get rounded up values in my computer.

For example, this is the webtool output:

#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1  100.00   92.79   65.29   63.98
     2: M97264.1     92.79  100.00   65.66   63.92
     3: AJ292230.1   65.29   65.66  100.00   66.88
     4: KM379144.1   63.98   63.92   66.88  100.00

And this is the output in my computer.

#
#
#  Percent Identity  Matrix - created by Clustal2.1
#
#

     1: JX310755.1     100      93      65      64
     2: M97264.1        93     100      66      64
     3: AJ292230.1      65      66     100      67
     4: KM379144.1      64      64      67     100

I have both installed clustalw using sudo apt install clustalw and by downloading the tarballs from http://www.clustal.org/download/current/clustalw-2.1.tar.gz, with exact results.

The commmand I'm using is clustalw2 alignment -tree -pim.

Any thoughts on why this could be? Thanks in advance.

Try clustal omega, instead of clustalw because it is really old. It should report to 2 dp. Variations do occur depending on how indels are treated, but this is just rounding.

You could try

alignment -tree -pim -score

The idea is to convert to absolute rather than percentage and might overcome the rounding error


How to get transcriptome FASTA file for viruses for Kallisto pseudo-alignment? I guess the title is self-explanatory. I'd like to find the reference transcriptome files for a few human viruses. My RNA-seq samples are from human tumor cells that were infected with viruses. The main part of my analysis is quantifying human gene expression from these samples. However, I also would like to find what proportion of my raw data is originated from virus RNA.
So, I'd like to learn how to find the reference transcriptome data for human viruses. Moreover, I appreciate it if you could give me some insight into your workflow for approaching this problem. Finally, if there is any specific command line, R, or python tool for this specific task, I'd like to know about it.

You might be interested in the viral sequence data from CATCH. There are various different viral subsets available.
You could also try the NCBI virus website, filtering on human hosts. There's a "Download" link on this page that allows you to download the full sequence set:
https://www.ncbi.nlm.nih.gov/labs/virus/vssi/#/virus?VirusLineage_ss=Viruses,%20taxid:10239&HostLineage_ss=humans,%20taxid:9605&SeqType_s=Nucleotide
You've asked a very specific question, so I'm giving you a very specific answer. It seems that what you actually want to do is different from the self-explanatory question.
For the purpose of creating an index for Kallisto or Salmon, it should be okay to use the full dataset. Those programs are designed to handle eukaryotic transcriptomes as input, which typically include many gene copies and isoforms, so transcript / genome duplication shouldn't be much of an issue.
Without knowing more about your particular application, it's difficult to provide additional advice. From what you've written, i.e. "what proportion of my raw data is originated from virus RNA", a hit to any viral sequence would be sufficient. For example, if you want any more specific information about viral proportion / species, then you should be using a metagenomic tool like Kraken2, rather than Kallisto / Salmon.


Output from vcftools missingness I'm new to data filtering on vcf data and vcftools.
I performed variant calling on my dataset, CHR22, homo sapiens. I'd like to remove sites that are missing in more than 5% of individuals.

vcftools --missing-site --vcf updated_ids68.vcf

This gives me a file with out.lmiss. That is missingness per locus.
I had 68 individuals in my input vcf file.

CHR     POS     N_DATA  N_GENOTYPE_FILTERED     N_MISS  F_MISS
22      16848278        106     0       30      0.283019
22      16848492        68      0       68      1
22      16849180        69      0       67      0.971014
22      16849229        68      0       68      1
22      16849376        133     0       3       0.0225564
22      16849476        132     0       6       0.0454545
22      16851734        126     0       36      0.285714
22      16852588        123     0       13      0.105691

I'm unable to understand the output here. What does third column - N_DATA tell?
133, 126 doesn't tell that this is number of individuals if my input had 68 individuals only.

VCFtools (0.1.15)
fileformat=VCFv4.1

Humans are diploid, so you can expect to see up to 2*N (2*68=136) alleles, so N_DATA is the number of observances of that allele.


When should I apply a multiple comparison adjustment? I am performing a RQL analysis ("http://www.esapubs.org/archive/ecol/E095/002/suppl-1.pdf").
This analysis meassures the relationship between different traits of three different matrix showing if there is a positive or negative relationship. This RQL analysis make a fourth-corner method.

Well, this method have the option to  make a p-adjust method to solve multiple comparison problem. However, my variables are independent and I belive that I should not apply to my dataset. Am I right? Or always that two or more matrix are compared I should apply a multiple comparison adjustment?

Thanks in advance.

The p-value adjustments corrects the inflated number of significant results when performing many significance tests.
If you have a test with 5 % false positive rate and but you run it several times the chances of getting a false positive results in any of these tests is far greater than 5 %. For example repeating the test with 5 % false positive rate 20 times means that there is a 60 % chance of at least one false positive result.

You can see this in the tutorial you have linked to in your question, where the number of significant differences goes down from 26 without correction to 18 with correction.

Therefore, if you are performing multiple comparisons you should do a correction. The problem of inflated p-values will occur when you are performing many tests even if the tested variables are independent.


You can read about it on the website, the last column is the Log likelihood scores (LLS) between the genes in the first two columns. An LLS of 0 means not better than random link.

On the website they show papers with fuller details about the statistics.


Is it possible to predict protein-ligand binding kinetics using machine learning? I would like to work on a project that involves the prediction of protein-ligand binding kinetics. What might be the feature that is relevant for the prediction?

This is a comment to redirect to other answers, but go too long.

Docking is the in silico prediction of where a ligand binds in a protein.
There are several Q&As about docking here that are relevant, so please check out:

machine learning, docking and the factors that affect it: Generate ligands candidates based on protein shape —also worth a read the "learning from the ligand" paper paper, which discusses why some ligand only parameters tell a lot about binding.
correlation between docking scores and enzyme kinetics: Purpose of correlating docking score with catalytic efficiency(Kcat/km) —discusses a paper https://www.nature.com/articles/ncomms12965 wherein a linear regression is done between docking scores and kcat/KM


Shannon's Entropy and DNA I have been working with Shannon's entropy (SE) and DNA lately, and referring to the formula and concept of SE, I just wondered whether one should in the case of DNA use each nucleotide as having a probability to show up as 0,25 or if it was better to use the probability based on the frequency in the genome (let's say the human genome).
In the case of the latter being the better, is it interesting to have those frequencies at the scale of the genome, the chromosome, the sub segment ?
I'm interested in calculating the entropy relative to a short window (the length of an Illumina read, about 125-150bases), using the proportion of each base on a given sequence length. The idea is "how much information is given" in this sequence, and by information I mean how much different bases do I have.
It might that I'm missing some basic mathematics knowledge here, don't hesitate to point it out so I can understand it better
cheers,

I think that genome-wide probabilities (i.e. not 0.25 for each base) should be fine. For a given segment of DNA, the chance that the sequence doesn't appear anywhere else in the genome is fairly high.
There is lots of patchy non-randomness (e.g. at centromeres, STRs, and other tandem repeat regions), but there's not much point in modelling that because it's not random.


How can I transform a mapped BAM file into an unmapped BAM file? In order to use MergeBamAlignment (Picard), I need an unmapped BAM file and a mapped BAM file. I have two mapped BAM files:

one with reads mapped to the reference I want but without metadata such as read group
a second one with reads mapped to the incorrect reference but with metadata

When I try to use MergeBamAlignment with these two files, I get errors, such as "Second read from pair not found in unmapped bam". So I am trying to unmap the reads from my second BAM file. How can I do that? Is there a tool for it? Otherwise, I'll have to manually change the bitwise FLAGs, reference names, CIGARS, etc.

Looks like you need Picard's RevertSam operation (described in context here, manual here).

java -jar picard.jar RevertSam \
     I=input.bam \
     O=reverted.bam


Why my bim file doesn't match to my ped file as the Plink documentation suggests? Plink documentation about .bim says the columns 5 and 6 are the Allele 1 and Allele 2 respectively.
The documentation about the .ped in turn, says "...The first six fields are the same as those in a .fam file. The seventh and eighth fields are allele calls for the first variant in the .map file ('0' = no call); the 9th and 10th are allele calls for the second variant; and so on.."

I've simulated a dataset as per:

$ plink1.9 --simulate ds1.sim acgt --out ds1
...
$ head ds1.bim
1       null_0  0       1       C       T
1       null_1  0       2       A       C
1       null_2  0       3       C       A
1       null_3  0       4       T       A
1       null_4  0       5       G       C
1       null_5  0       6       T       C
1       null_6  0       7       T       C
1       null_7  0       8       G       C
1       null_8  0       9       C       G
1       null_9  0       10      T       A

Than I converted it to plain PED and MAP:

$ plink1.9 --bfile ds1 --recode --out ds1
...
$ head ds1.ped > ds1.head.ped
$ vim ds1.head.ped
per0 per0 0 0 2 2 T T A C A A T A C C T C T C C C....

Given the data above, why my .bim file shows "C and T" for the first variant while my .ped file shows "T and T"?
Apparently they don't match as per described in the documentation (or more likely, I'm misinterpreting it)...

Also, how do those 2 files related to each other? E.g.: How do I know that my Sample 1 has the SNP's null_3, null_4 or null_5, etc?

I need to clarify those doubt because I want to fit a LightGBM model with the Genomic Data (can be plink or vcf). Not sure yet how to do. I was trying to convert the plink files to a CSV-style somehow, to be able to use in my model. Not sure if this is the way to go though...

Thanks!

Responded in the plink2-users group (https://groups.google.com/forum/#!topic/plink2-users/zgJxdXxvdLo ).


What is the meaning of the "*" character in bwa fastmap's output? I am mapping kmers back to a few bacterial genomes using bwa fastmap:

bwa fastmap -l 9 ref.fasta kmers.fasta > out.fastmap
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[main] Version: 0.7.17-r1188
[main] CMD: bwa fastmap -l 9 ref.fasta kmers.fasta
[main] Real time: 0.168 sec; CPU: 0.157 sec

The output looks like this:

SQ      AAAAGTAGTAAGCAGGAAGACAACACGGTTG 31
EM      0       11      1       contig1:+2041368
EM      3       13      2       contig1:+1875252        contig1:-2474744
EM      4       15      1       contig1:+3779779
EM      5       16      1       contig1:-3253348
EM      6       17      2       contig1:+938710 contig1:+1066682
EM      7       20      1       contig1:-4912797
EM      10      21      4       contig1:+4803473        contig1:-2830252        contig1:+4495150        contig1:-4907283
EM      11      23      1       contig1:-3148132
EM      13      24      1       contig1:-4651172
EM      14      25      3       contig1:-2142223        contig1:+4994474        contig1:+873066
EM      16      26      4       contig1:+156775 contig1:+27749  contig1:+2207492        contig1:-1340811
EM      17      29      1       contig1:+3523989
EM      19      30      24      *

EM      20      31      3       contig1:-2533354        contig1:-208660 contig1:-1080177
//

My question is: what is the meaning of the * char in the row before last? Why is bwa fastmap not reporting all the hits as in the other lines? I could not find a page that explains this output anywhere.

It means the read is unmapped. No origin in the reference could be found that was sufficiently similar to the read to call it a proper alignment.

In the SAM specification, * is always used if the information for that field is not available. Therefore, if no mapping information is available (=unmapped) one sets *.


10x Genomics Chromium single-cell RNA-seq data analysis options? Provide an overview of 10x data analysis packages.

10x provides Cell Ranger which prepares a count matrix from the bcl sequencer output files and other files (see bottom of page https://support.10xgenomics.com/docs/license for the programs it uses).

What can we do with the output files?

Data preparation

Cell Ranger uses the Illumina sequencing output (.bcl) files

Make fastq files:
cellranger mkfastq ==> .fastq
Prepare count matrix: cellranger count ==> matrix.mtx, web_summary.html, cloupe.cloupe
Optional: combine multiple matrix.mtx files (libraries): cellranger aggr

Data analysis

Loupe Cell Browser visualization of cloupe.cloupe files

Count table matrix.mtx analysis options:

Python
R Cell Ranger R Kit: cellrangerRkit::load_cellranger_matrix() ==> ExpressionSet
R Scater: scater::read10XResults() ==> SCESet object
R Seurat: Seurat::Read10X() ==> Seurat object


Using R to import specific records from a PDB file I'm using R and I'm trying to import specific records from a PDB file (e.g. HETNAM, HETSYN). Unfortunately, it seems the read.pdb function of either biomaRt or Rpdb aren't capable of it.

Any ideas?

Thanks

If your file is in PDB format as the one shown here in .pdb format then you can simply extract your desired fields using sed which is a programming language used to parse and transform text. As for this structure i.e 6C4A in .pdb format HETNAM and HETSYN can be extracted using sed as follows:

sed  -n -e '/HETNAM/ p' -e '/HETSYN/ p'  6C4A.pdb

HETNAM     EJA S-[(1Z)-2-CARBOXY-N-HYDROXYETHANIMIDOYL]-L-CYSTEINE
HETNAM     PYR PYRUVIC ACID
HETNAM      MG MAGNESIUM ION
HETNAM     3NP 3-NITROPROPANOIC ACID
HETNAM     PEG DI(HYDROXYETHYL)ETHER
HETNAM     EPE 4-(2-HYDROXYETHYL)-1-PIPERAZINE ETHANESULFONIC ACID
HETSYN     EPE HEPES


Using the t-SNE algorithm on microarray data + an error bonus I'm trying to use the t-SNE algorithm on some microarrays data. More specifically my data frame has 18600 columns with genes (features) and 72 rows with conditions with replicates ( 10xWt , 10xTg , etc ). The expression values are in log2 scale.

Here is the code that I'm trying to run.

# t-SNE implementation
library(Rtsne)

set.seed(1)

for(i in 1:15){

  tsne = Rtsne(data.T[,-18601], dims = 2, perplexity=i, verbose=TRUE, max_iter = 1000, pca=T)

  colors = rainbow(length(unique(data.T$classes)))
  names(colors) = unique(data.T$classes)
  plot(tsne$Y, t='n', main="tsne")
  text(tsne$Y, labels=data.T$classes, col=colors[data.T$classes])

  readline(prompt="Press [enter] to continue")
}

Please note that I'm not counting the column 18601 because this colums contains the labels/classes for each condition.

The think here is that when I execute this script, R returns me this error:

  Error: protect(): protection stack overflow

Should I change the --max-pp-size or it's a bug in Rtsne package?

Also I was wondering if it is more meaningful to run the tSNE algorithm using not the log2 values of the expression level but the log fold change values in respect to the Wt (wild type) condition. I'm asking because I couldn't find a such other implementation of the tSNE on microarray data.

For the configuration of the Rtsne function I read this article

Any other suggestion on the implementation is welcomed.

Converting your data.frame to a matrix (and then removing the data.frame) will often free up enough memory that you won't run into this. Note that a matrix is more memory efficient than a data.frame and you're requiring Rtsne() to hold both in memory at the same time (many math-centric functions will end up converting things to a matrix at some point for efficiency).

For what it's worth, it's never been entirely clear to me what the interaction is between a data.frame and the pointer protection stack, but it's often the case that this solves this sort of error.


Protein PDB/CIF text file reader that can fold a same section into one line Is there a software that can recognise the protein PDB/CIF files, and fold multiple lines of a same section into one line, so that it is easier to handle? I need to understand the format in CIF so as to create my own CIF file.

For example, when Notepad++ is used to read a .py file, it can fold the same section into one line.

So for a PDB/CIF file, e.g. all the atom lines for a same chain can be folded into one line, and I can expand it if needed.

Like this picture, assuming "This PC" is the whole protein, while the other objects are the different chains.

Programs that work with PDB and mmCIF files may present the structure as a tree, for example:

(this is a Find dialog from ccp4MG)

but I'd be surprised if a text editor had such a feature. Normally, you don't edit PDB files in a text editor.


How I deal with this kind of gene expression comparision I have RNA-seq from two sequencing batches; Lab technician says that he has run the RNA expression quantification two times in bathes 1 and 2 for example tumor 1 in batch 1 and tumor 1 in batch 2 , normal 2 in batch1 and normal 2 in batch 2. PCA shows however shows two experimental runs are too close to each other and can be merged

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

Now my question is can I compare only a tumour sample with its own matched normal sample in one patient to get differentially expressed genes for this patient by using experimental runs (technical replicates )? Design would be Tumour1_batch1 Tumour1_batch2 versus Normal1_batch1 Normal1_bacth2

> head(mycols)
              Condition
Normal1_batch1  N
Normal1_batch2  N
Tumour1_batch1  T
Tumour1_batch2  T

No. If the only difference between Tumour_batch1 and Tumour_batch2 was the library prep/sequencing run, then all this will tell you is if the particular sample of RNA taken from the tumour is different from the particular sample of RNA taken from normal. It almost certainly is.

Indeed if you did technical replicates of two RNA samples, both taken from normal, you would almost certainly find differences.

What you want to find is the differences in RNA level between tumour in general and normal in general. To do that you need to make an estimate of the mean and variance of RNA levels in the tumour and in the normal. To do that you need multiple samples from each.

If you need differences within a single patient, it would be more honest just to compute the fold changes without any implication of statistical certainty. You would then need to think about alternate schemes for measuring the confidence in inferences drawn from the data when you consider the actual biological questions.


"Sequence Duplication Levels" module still fails after pre-processing Illumina data I want to ask about why the sequence duplication levels are high after I trimmed by using Trimmomatic? I am using the following Trimmomatic operations: HEADCROP = 19 TRAILING = 20 MINLEN = 66.

How can i solve this problem? Thank You.

FastQC assumes that all samples are for whole genome sequencing and will flag them as failed if they differ too much from that assumption. This will, for example, cause essentially all RNA-seq, ChIP-seq, and ATAC-seq samples to fail in one module or another. This is not any cause for concern and is completely expected. Primarily concern yourself with whether all of your samples are similar in their metrics.


How to retrieve the DNA sequence from a particular species and region through command line? I am new to working with the NCBI database so I am not sure how trivial this question is. I wanted to get DNA sequence of a specific version of the Xenopus tropicalis genome of a specific region (ie 1-100bp) without doing it manually each time (going on genome browser). I know there is a faster way and efficient way to do this through the command line but I found the NCBI tutorial hard to follow.

Samtools
Samtools is a toolkit that includes a binary called samtoools, which can be used to query a FASTA file — to extract a subsequence based on some coordinate range you specify.
If you do not have samtools installed, you need to add it to your computer. For OS X, you could use Homebrew, for example:
$ brew install samtools

For Ubuntu, you could use apt-get:
$ sudo apt-get install samtools

Reference assembly
The next step is to download a reference assembly for your genome and genomic build-of-interest.
There is a research group dedicated to making Xenopus datasets available. If you were interested in v10 of the Xentr genome, for instance, you could do the following or similar to download it:
$ curl --output "XENTR_10.0_genome.fasta.gz" "http://ftp.xenbase.org/pub/Genomics/JGI/Xentr10.0/XENTR_10.0_genome.fasta.gz"

You could also dig through the NCBI site to get a similar link. It's the same idea — just change the web address, accordingly.
Indexing
Indexing makes it possible to extract an arbitrary subsequence of interest from this FASTA file. You need to index, but you only need to do it once:
$ samtools faidx XENTR_10.0_genome.fasta.gz

This creates two new files: XENTR_10.0_genome.fasta.gz.fai and XENTR_10.0_genome.fasta.gz.gzi.
Note: You need to keep these two files in the same directory as the original, compressed FASTA file XENTR_10.0_genome.fasta.gz.
Querying
Now you are ready to query, or extract a subsequence of interest. If you wanted the first hundred bases of the first chromosome, for instance:
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100
>Chr1:1-100
cctcagcgtcttcttctttctgctcccactgcattagcatagggagagagggccacagag
gctaggtagcattgccggctagcaatcttcagcgtcctct

The search term used here is Chr1:1-100. Replace Chr1 with the chromosome of interest and 1 and 100 with the start and stop coordinates, respectively, of interest. Coordinates used with samtools queries use a one-based indexing scheme.
Note: This uses the chromosome naming scheme in the XENTR_10.0_genome.fasta.gz file (Chr1, etc.). Other research groups (NCBI) might use different chromosome names. If you get your FASTA from somewhere else, investigate the FASTA file to determine the naming scheme being used.
If you want to redirect output to a file, use the output redirect operator (>):
$ samtools faidx XENTR_10.0_genome.fasta.gz Chr1:1-100 > mySubsequence.fa


Large dataset normalization for PCA I need to normalize a large (400Mb) dataset for doing PCA analysis. I want to use scran for doing that:

biocLite("SingleCellExperiment")
biocLite("scran")
library(SingleCellExperiment)
library(scran)
list_of_sce <- list()
# Looping though the UMI_count 'split_factor' columns at a time
split_factor = 500
for(i in seq(1,ncols, split_factor)) {
  num_loop = floor(i / split_factor) + 1
  idx = ncols
  if (i + split_factor < ncols) {
  idx = i + split_factor
  }
  sce <- SingleCellExperiment(list(counts=UMI_count[, i : idx]))
  # Normalization of the dataset containing
  # heterogenous cell data (different cell types)
  clusters <- quickCluster(sce) ## <- ** error appears here **
  sce <- computeSumFactors(sce, cluster=clusters)
  list_of_sce[[num_loop]] <- sce
}
cbind(list_of_sce)

However, when I try to run this it produces an error at the quickCluster step:

Error: cannot allocate vector of size 23.7 Gb

It feels like cbind at the end is not the right way to go, and I need to somehow combine sce and re-normalize them again somehow. The question is how to do that in a correct way: normalizing large dataset for PCA analysis when it can not be fully loaded into memory?

I normalized a high throughput dataset for a school project using DESeq library using the script bellow. The code is based on a lesson I had. My goal was determine the over expressed genes, but the normalization step should be the same.

#####################################################################
### DATA PRE-PROCESSING: Normalising gene expression distributions
#####################################################################
# Genes that are not expressed at a biologically meaningful level
# in any condition should be discarded to reduce the subset of genes
# to those that are of interest, and to reduce the number of tests carried out
# downstream when looking at differential expression.
# Using a nominal CPM value of 1 (which is equivalent to a log-CPM value of 0
# a gene is deemed to be expressed in a given sample if its transformed count
# is above this threshold, and unexpressed otherwise.
# Genes must be expressed in at least one group (or roughly any three samples
# across the entire experiment) to be kept for downstream analysis.

library("DESeq")

keep.exprs <- rowSums(cpm>1)>=3 # is filtering is ok
DGEfiltered <- DGEobject[keep.exprs,, keep.lib.sizes = FALSE]

## Normalisation generally required to ensure that expression distributions
## of each sample are similar across the entire experiment.
## Any plot showing the per sample expression distributions, such as a density or boxplot,
## is useful in determining whether any samples are dissimilar to others.
## Normalisation by the method of trimmed mean of M-values (TMM)
## is performed using the calcNormFactors function in edgeR.
# cpmF = filtered subset # log.cpmF = log cpm filtered subset
# Copy unnormalized, filtered data to show in boxplot
DGEf.notnorm <- DGEfiltered # normaliseerde
DGEf.norm <- calcNormFactors(DGEfiltered, method = "TMM")
DGEf.norm$samples$norm.factors

# --> scaling factors relatively close to 1 --> mild TMM normalisation

## Boxplots of log-CPM values showing expression distributions
## for unnormalised data (A) and normalised data (B) for each sample
par(mfrow = c(1,1))
par(mar = c(7,4,2,1))
# Boxplot unnormalised data
DGEf.notnorm$samples$norm.factors
# [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1
log.cpmF.notnorm <- cpm(DGEf.notnorm, log = TRUE)
boxplot(log.cpmF.notnorm, las = 2, col = sample.col, cex = 0.9,
        main = "A. Unnormalised data", ylab = "Log-cpm")
# Boxplot normalised data
DGEf.norm$samples$norm.factors

log.cpmF.norm <- cpm(DGEf.norm, log = TRUE)
boxplot(log.cpmF.norm, las = 2, col = sample.col, cex = 0.9,
        main = "B. Normalised data", ylab = "Log-cpm")
# --> only when scaling factors are not close to 1 you will see a big difference


Why is ribosomal RNA difficult to remove even with Poly(A) selection? In this answer (actually in a comment), it is stated that:

  As you've noticed from your own analysis, the ribosomal genes have quite variable expression across cells. They're expressed everywhere, and quite difficult to completely remove from a sample (even with polyA selection), with success of removal depending on things like the amount of mechanical damage that the RNA has been exposed to.

Which are the most relevant issues that make Poly(A) selection fail? Looking online, Poly(A) is suggested often but I didn't find an explanation for the different ways in which it can fail.

We've found ribosomal RNA to be less of a problem with sequencing that depends on polyA, which suggests the issue might be in the library preparation, rather than the selection.

Many polyA RNA library preparation methods involve amplification, rather than selection, which means that existing transcripts that are present in very high abundance (such as rRNA) will still show up in the mixture.

But polyA-tailed RNA isn't the only RNA present in a sample. For a comprehensive look at RNA, the whole-RNA sample preparation methods can't use any form of selection to identify target molecules. In this case, it's typically recommended to use a form of ribosomal RNA removal, where a probe to ribosomal RNA is attached to magnetic beads (or a column, or similar), and used to fish out out ribosomal RNA, leaving the remainder behind. Unfortunately those probes don't work perfectly, either due to the ribosomal RNA being damaged at the probe site, or the probe not matching the ribosomal RNA molecule, or the RNA just not getting into the right physical state to bind with the probe. All these things mean incomplete removal of ribosomal RNA, ending up with a situation where the highly-abundant ribosomal RNA still pokes through in the reads.


Splitting a GenBank file into smaller files I have a genbank file that I have been trying to run in HMMER using the following command:
hmmscan --domtblout NC_000962.3.final.gbk.out -E 0.000001 prot_profile.hmm NC_000962.3.final.gbk

However, since the file is too big for HMMER, I get the following error:
Fatal exception (source file p7_pipeline.c, line 697):
Target sequence length > 100K, over comparison pipeline limit.
(Did you mean to use nhmmer/nhmmscan?)
Aborted (core dumped)

I need to split this file into smaller files for the HMMER run. The problem is, I would need the ORIGIN field (containing the nucleotide sequence) at the bottom of the file split into relevant chunks into the subfile as well. This is where I am getting stuck.
I have already looked at GenBank Parser and seqretsplit to no avail.
I would greatly appreciate some help with this. Thank you.

You can use NCBI EntrezDirect to download data in GenBank format for a specific region of the sequence as follows:
efetch -db nuccore -id NC_000962.3 -format gb -seq_start 1 -seq_stop 99000 > seq.gb


A box plot of qualitative variables I have 2 groups of patients; responders to chemotherapy and non-responders

I have calculated cancer cell fraction (CCF) of a set of driver genes (by variant allele frequency) for each group and I have

> head(dat)
    Response       CCF
1 Responders 1.0000000
2 Responders 0.5413323
3 Responders 1.0000000
4 Responders 1.0000000
5 Responders 1.0000000
6 Responders 1.0000000

unique(dat$Response)

[1] Responders    Nonresponders
Levels: Nonresponders Responders
>

If CCF > 0.95 , the mutation is clonal otherwise sub clonal

I want to show how many clone and subclone are in these two groups by box plot or something similar like below

I have tried this which was nonsense

ggplot(dat, aes(x=Response, y=CCF)) +
    geom_boxplot()

Can you help me?

I'm not sure that boxplot will be the more appropriate representation, as you will end with two numbers (count of Clonal and SubClonal) per groups of patients.

One solution will be to first create a new categorical variable based on CCF values using for example ifelse statement for example by writing:

df$Clonal <- ifelse(df$CCF > 0.95,"Clone","SubClonal")

Then you can get the count of each Clone and subclone for each responders and non responders by using table:

DF <- as.data.frame(table(df[,c("Response","Clonal")]))

You can finally convert these new values as percentage by grouping according to the Response. For example, using dplyr, you can do something like that:

library(dplyr)

DF %>% group_by(Response) %>% mutate(Freq_percent = Freq / sum(Freq))

And finally, you can plot it as a barchart or single point in ggplot2.
Hope it helps you to figure it out how to deal with your data.


Excess average estimated identical-by-descent in genotype data (Cross-post with Biostars: https://www.biostars.org/p/470788/)
We have some genotype data that we are putting through quality control in PLINK 1.9. As part of this QC, we have limited the data to subjects of self-reported and genomic European ancestries, and have run basic QC steps, removing call rates below 95% for individuals and variants, and removing variants excessively out of Hardy-Weinberg equilibrium (< 10^-10).
Pruning for LD and limiting to the autosomes, we have calculated estimated IBD between our participants and the estimated the average pairwise IBD for each participant. We have found that an unusual number of our participants (233/13049) show a high mean pairwise IBD with others (>3SD from the cohort mean, > 0.028 in this instance, maximum 0.075). Examining examples of the individuals in question shows that this is driven by low-level relatedness (~0.03-0.1) with many other individuals, rather than a small subset of relatives in the data. In addition, the relatedness in each case stems entirely from IBD1 (one shared variant) and not at all from IBD2 (two shared variants). Individual average IBD correlates negatively (-0.69) with the genome-wide inbreeding coefficient.
I have one idea about what could be causing this (partial cross-contamination of the samples), but I am interested to know if others have other ideas what would be the cause?

A clear negative association within an inbreeding coefficient suggests convergent evolution driven by shared selective pressure. There is an outside chance of purifying selection, or else balancing selection, but selection is the key word. The phenotype of the locus needs investigating, i.e. outside nearly neutral theory what is constitutes positive selection?
Allopatry caused by genetic isolation and genetic drift is your null hypothesis. However Darwinian positive selection is the alternative hypothesis.
I would be wary against suggesting contamination unless you've a clear suspicion this is occuring. I do accept that in human genetics its often difficult to avoid, particularly in ancient DNA type work. In this case you would need to try and control for it. The argument I would put against contamination is that 'genetic relatedness' is localised within the genome, whereas contamination would be a genome-wide phenomena.
Convergence is common within for example HLA genotypes (T-cell) loci between global populations.
You could perform formal tests for postive selection either within pop gen (the locus will associate with significant homozygosity) or phylogeny (point mutation).


Determine if an organism is microbic or unicellular I am doing some metagenomic analysis of samples which can only contain microbes due to the experimental setup.
For this, it would be useful to be able to distinguish whether a specific organism is a microbe or not, as it allows:

to detect potentially false findings,
to avoid an unnecessary inflation of the metagenomic database.

For animals and plants, this is straightforward, as they are never microbes.
However, for fungi (and some other clades) it gets complicated.
At the very least, I failed to find a simple list of microbial or non-microbial clades.

Ideally, I would like to have a function that maps a taxonomic ID to a boolean stating whether this organism is microbial. I acknowledge that this is utopic since microbe is not that well defined, there are facultatively multicellular organisms, and for some cases it may simply not be known. I am therefore happy with any reasonable approximation such as unicellularity. In particular, given my application, a small chance of falsely identifying an organism as microbial is not an issue.

I am currently working with the NCBI database accessed via the ETE 3 module for Python, but I am open to use other databases and tools.

I would start by saying that your statement that no animals/plants are microbes is false. There are indeed microbial algae (plants) and microbial animals.

So it is indeed difficult. (I'll briefly note that unless your protocol specifically gets rid of environmental DNA, you will probably have some contaminating non-microbial DNA no matter what.)

However, we can use a different version of your proposition to make a simple (over-simple, most would argue) approximation: all prokaryotic organisms, e.g. bacteria and archaea, are microbial, and all eukaryotic organisms are not. This is not hard to do if you are already using NCBI taxonomy, as you can use e.g. this mapping. It has the advantage that all bacteria and archaea are definitely microbial.

Of course it's quite wrong, as you already know. In principle you could set up a complex decision tree sorting this out, though I imagine that's the thing you wish to avoid.

In terms of other genome databases, you might consider looking into GOLD, which has extensive metadata in a slightly less arcane setup than NCBI, which might be helpful. But there is a high rate of missingness there, and also there is still not necessarily a great "microbe"/"not microbe" mapping.

There are also existing databases of microorganisms, which may or may not satisfy your requirements; I have little familiarity with most of those so caveat emptor.

I believe (hope) that NCBI taxid is a good way to do lookups in all such databases, but it may not be true.


Question about research in bioinformatics for Computer Science student I am a computer science student, who was interested in taking research in the field of bioinformatics.

The first research idea I had was to make a model that could predict the possibility of type 2 diabetes by using machine learning such as convolutional neural networks or recurrent neural networks by training the model with fasta data of protein that causes diabetes.

The research will focus more on the architecture of the machine learning itself.

Does this research possible to work on ? If yes, how and where can I get the data?
All my lecturers have no background in bioinformatic, so I can't ask any of them.

Edit:
I m sorry, not protein but gene like CAPN10, CDKAL1, CDKN2A, HHEX, HNF4A, IGF2BP2, KCNJ11, PPARG, SLC30A8, and TCF7L2 that shown have risk to be associated in Type 2 Diabetes (http://diabetes.diabetesjournals.org/content/56/12/3033).

I think , this kind of research never stop until find way for cure .
I have some research in heart disease,and write this paperECG, you can have some research in improve learning method to find fast and more accurate diabetes , there are lot of database your can search in google and more famous is UCI   .


Get list of urls of GSM data set of a GSE set I have this GSE dataset ( GSE104279 ) (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE104279).
I want to make a table with set IDs and ftp urls to use it as a table in galaxy.org
I know that we can use ENA to get specific arrangement:
https://www.ebi.ac.uk/ena/browser/
I tried to get:
SampleID    Group   URL

so I used :
https://www.ebi.ac.uk/ena/browser/view/PRJNA412223
But nothing is showing.
Is there away to get these urls in the arrangement above?

You can use Entrez Direct for this as follows:
esearch -db gds -query 'GSE104279' \
  | esummary \
  | xtract -pattern DocumentSummary \
    -if 'entryType' -equals 'GSM' \
    -def 'NA' -element Accession title summary FTPLink

This will return a table with data similar to this:
GSM2580330      18_Z13_2_d5_Zika        cortical organoids_Zika_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580330/
GSM2580329      17_Z13_2_d5_Control     cortical organoids_mock_5d      ftp://ftp.ncbi.nlm.nih.gov/geo/samples/GSM2580nnn/GSM2580329/

To download sequence reads, you should follow links to SRA. Using Entrez Direct you can do this as follows:
esearch -db gds -query 'GSE104279' \
  | elink -target sra \
  | efetch -format runinfo

This will return a comma-delimited table containing SRA identifiers and an FTP path to the SRA data. These FTP paths won't lead you to FASTQ files though. You can pass the SRA run identifiers of the format SRR### (or ERR### or DRR###) to fastq-dump or fasterq-dump tools from the SRAToolkit to download data in FASTQ format.


Biohackers Netflix - DNA to binary and video I'm not sponsored or anything, just interested in their challenge to decipher their DNA code.
They encoded their first episode of "Biohackers" video/binary file to DNA code and said if we could decode it we can watch it (without Netflix).
Here's their page: https://biohackersnetflix.com with description and download for the DNA sequence file.
(Don't know if it's just in German or you can translate it. If questions regarding this page, ask me.)
The file is ~550MB small and contains 3.882.771 lines (not in fasta format). Every line has a length of 147 characters including primers at both ends (Illumina?).
Here are the first 5 lines:
ACACGACGCTCTTCCGATCTCTCCCAGGGACAAAGGTTCTGCATTTGCAGCAAGACTCCTGTAGTGCTGCAGATTCTCTGGTTGGATAGTACGGCGTACATTTCTGTATTGTAGCACCATGGGGTAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTAAGGCTTCGTAACAGATATTCTATATCGTCACATTGGTCTGAAGGAAGTCGCCTATAATCGCTCCTCTGTTTTTTAAAACTGCTATGGACCCGCTGTTCGGTGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTCATGGTATAAGTGTTAAGGGTAATAACCACCTACCCCCCTCATTGCTCGTTTTTCCTGGAACCTTAACATTCGCAATAGCTAGCTGTTTCCTAGTAGAACCAAGGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTAGGATGTAGTCACAGGTCATTGTCATTAACTCAACCGAGGACATAACACTAAGTCCCACTAGGCCTGGATTCTCTAACGCGGTCTCTCTATTGGGGGAAGGGGTGAGATCGGAAGAGCACACGTCT
ACACGACGCTCTTCCGATCTTCTGGTAAGGCGGGTTGATATCAGTCACCTCCCTTTGAGCTAAAATACGATGGCGATTTAGTGTGAAACTAATAATGCTTGTCATACCAGCAGTACCGGATCGGGAGATCGGAAGAGCACACGTCT

I trimmed all the primers and tried to decode {A, C, G, T} considering every permutation {00, 01, 10, 11} as the obvious(?) decryption method (4! = 24 possible decodings) using python.
Then I hoped to get 1 of these 24 files loaded into VLC media player or something to be played, but it didn't work and every file seemed to be broken in the same way. I think I'm missing something here.
Can I assume that a text file containing only 0's and 1's should be playable in VLC if the DNA code is correctly decrypted?
(If I am wrong here, please tell or move me.)
//Edit: I converted all 24 files to ASCII to see if there's some kind of "video-like header". (All videos have some sort of description in their first lines if opened in text editor?) But there's just gibberish.
//Edit: I saw that every 84th sequence position has a "T", which is kind of weird.
So I tried to run my script again with these T's removed, but still no solution.
//Edit: I searched for "AVI", "264", "codec" and some other strings in every video file I created and hexdumped. Nothing found.
For clarification: I translated the DNA into every 24 binary and then into their ASCII representation following the 19 upvotes answer: https://stackoverflow.com/questions/7290943/write-a-string-of-1s-and-0s-to-a-binary-file. The 104 bases / 208 Bits (removed repetetive "T" and primers) are actually a multiple of 8 (respectively 26 Bytes) so I could be on the right way (even if not 32 Bytes?).
De novo Assembly didn't work and I found no obvious ORF "genes" representing some kind of URL to the video or something which was a neat idea considering the video file would be only ~150MB. (See comments.)

I was quite curious when I saw that as well. I've spent more time than I'd care to admit trying the sort of things you have. I've gotten it decoded now, but I can't really claim any sort of victory, as the problem was pre-solved.
After struggling through some of the same experiments you did, I decided to take a closer look at their explainer video here: https://youtu.be/DMYgjOHgHxc
First I tried to decipher the blackboard sketch around 1:30, but it was quite vague. A closer investigation into the speaker, however, yielded some luck. Googling Dr. Reinhard Heckel brought me to his website here: http://www.reinhardheckel.com/ which shows his most recent publication – a paper on encoding digital data in DNA.
The encoding is relatively complex (something nearly impossible to stumble upon accidentally, I'd reckon) but fascinating. The fragments are indexed and have two, layered error-correcting codes. Perhaps making our job harder (but serving the practical purpose of minimizing homopolymers and unwanted annealing) is the fact that the data is XORed with a pseudo-random noise that shuffles the data. The paper obviously has all of the details.
From my skim, however, I stumbled across this Github repo linked in the paper: https://github.com/reinhardh/dna_rs_coding
As it happens, the README was just updated in this last week, describing how to decode the episode from the file provided. If you have Docker, you just need to copy-paste some commands.
The product is indeed the full and final episode – coming in at a dainty 63.1MB. How did they get 40 minutes of video that compressed? Well, in short, 720x360 at 24 fps. It looks a little awful and is in German only, but is certainly a cool little easter egg. If nothing else, I have a cool paper to read.


how to filter a multi-VCF for sex-specific genotypes? I have a single VCF file of variants for a population that includes males & females and I would like to pull out all variants that are sex-specific. For example, sites where all males are 0/1 and all females are 0/0. Is there a way to do this using bcftools or vcftools? I could enumerate each individual name, or ask for all individuals named *.M and *.F.

I first tried running smoove on the male and female populations separately then comparing the two VCFs, but I'm working with a draft genome with loads of contigs. As only contigs containing variants are output in each VCF, and these are different between the two files, vcftools couldn't handle the comparison.

Someone here https://github.com/samtools/bcftools/issues/118 managed to get it working ike this:

$bcftools view -f PASS -i 'GT~"0/0.*0/1.*0/1.*0/1.*0/1"' combined.vcf.gz

but the same syntax gives me no variants (yet I can see some that meet my criteria in IGV).

Found a great tool with SnpSift!

From the examples: I want to keep samples where the genotype for the first sample is homozygous variant and the genotype for the second sample is reference:

cat variants.vcf | java -jar SnpSift.jar filter "isHom( GEN[0] ) & isVariant( GEN[0] ) & isRef( GEN[1] )" > filtered.vcf


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

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

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

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


Number of sing-copy genes identified by orthofinder is far less than BUSCO assessment? Dear Orthofinder users and BUSCO users,
Recently, I have assessed 30 genome annotataions using BUSCO with lineage dataset ascomycota_odb10, I got 1100+ complete and single-copy genes. But the orthofinder  only identified 300+ single-copy genes. Aneone else have encountered this similar problem?
Thanks,
Yanpeng

It's hard to say for sure without seeing the exact commands you used. However, Orthofinder identifies/defines single-copy orthologues as those genes present in a single copy in all species being compared in the analysis. If a BUSCO orthologue is present in more than one copy in one of the species being compared (could be multi-copy, could be fragmented etc.), then it won't be recovered as single copy by Orthofinder.


What is the R/Bioconductor preferred way of accessing TCGA data? It seems there are two packages providing access to TCGA data these days:

TCGAbiolinks (older)
GenomicDataCommons (newer)

Will GenomicDataCommons supersede TCGAbiolinks? Is there anything that I might miss if I decide to use only GenomicDataCommons?

There isn't a R/Bioconductor person or view, so there isn't a preferred way to access TCGA. Use what suits you.

Bioconductor is build by contributors which are encourage to not overlap in their contributions. Usually the packages have different options/features even if they allow to access to the same data.

In this case, the TCGAbiolinks package is specifically for TCGA dataset while GenomicDataCommons is for "access the NIH / NCI Genomic Data Commons RESTful service". So the functionality they bring to R and Bioconductor is different.

The fact that one project is finishing doesn't mean that the data or the packages build around that project are obsolete. Some functionality different between them is the functions to harmonize and prepare (GDCprepare) data from the TCGAbiolink package, or see a boxplot correlation and AAIC plot to define outliers TCGAanalyze_Preprocessing.


Why does Coronavirus genome contain poly A region at the 3' end in GenBank? If you look at the Coronavirus genome below, it contains poly-A region at the end. I know that the viral RNA contains a poly-A tail. But, the genome in GenBank is stored in cDNA format. That means the end should have poly-T instead of A. Can anybody help me to understand this, please ?
>MN988713.1 Wuhan seafood market pneumonia virus isolate 2019-nCoV/USA-IL1/2020, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
--------------------------------------------------------------------
GATCGAGTGTACAGTGAACAATGCTAGGGAGAGCTGCCTATATGGAAGAGCCCTAATGTGTAAAATTAAT
TTTAGTAGTGCTATCCCCATGTGATTTTAATAGCTTCTTAGGAGAATGACAAAAAAAAAAAA

I figured out that it has to do with second strand of cDNA. Usually genbank sequences are equivalent to 5’ to 3’ strand of cDNA. This strand is generated using second strand synthesis enzyme.


How to make a video of a gromacs .trj file with VMD I did gmx trjconv on the GROMACS trr file to generate the .xtc file, and then used gmx dump to generate the trj file. Then I tried to make a movie with VMD using the trj file. As a result, it does not work properly. The error in the terminal is as follows.
gromacsplugin) Cannot read header fromm '/Users/kondouakira/mount/Lab/1ns_sim/md_0_1_noPBC3.trj', no error

I use following command.
gmx trjconv -s md.tpr -f md.trr -o md_0.trr -pbc mol -ur compact -b 1 -force
gmx dump -f md_0.trr > md_0.trj

No need to use gmx dump, you can simply load a xtc or trr trajectory in VMD using:
vmd path/to/gro_file.gro other/path/to/trajectory_file.xtc

As always you can also use the GUI to first open your gro file, then right click on the molecule and select Load data into molecule before selecting your xtc or trr file.
After putting the representations as you see fit, if you want to export this as a movie (GIF, MP4, others...) you can go in Extensions -> Visualization -> Movie Maker. See Trajectories and Movie Making VMD official tutorial for more details.


How many entries in PDB have SS bonds? Is there a way to figure out what percentage of entries in PDB have SS bonds?
In addition, what other types of bonds in PDB may be worth looking into (not H bonds)?

The easiest way is to check annotations in mmCIF or PDB files.
Looking at my page with statistics of mmCIF tags, last updated in March:
_struct_conn_type.id        59137×metalc 42706×covale 36195×disulf 10725×hydrog
_struct_conn.conn_type_id   5350402×hydrog 1733277×metalc 559092×covale 226472×disulf

Which means 226472 disulfide bonds in 36195 PDB entries at the time when the PDB had 161470 entries.


How to represent random sequence elements in SBOL? In our lab, we work with synthetic biology components with partially random sequences (similar to work in directed evolution). So, for example, we have a plasmid design with several components, including one protein coding sequence that will be created from a wild-type sequence via mutagenic PCR.
At the design stage, we don't know the exact DNA sequence for any of the resulting variants of that protein coding sequence. We just know the wild-type sequence and the average number of mutations we'd like to apply.
How can we represent that design in SBOL (via the Component class)?
Then, when we implement the design, we end up with a single sample that contains plasmids with 10^5-10^6 (or more) different protein coding sequences. We then measure the complete DNA sequence for every plasmid variant, but I don't want to end up with >10^5 different SBOL definitions to describe the sample.
How can we best represent that sample in SBOL (via the Implementation class)?
Finally, we also end up measuring the performance of every plasmid variant, and we identify the specific protein coding sequences for useful and/or interesting variants. We then design and build specific, clonal plasmids with those useful/interesting protein coding sequences.
Representing those clonal plasmids as SBOL Components seems pretty straight-forward, but is there a good way to capture the relationships with the original random-library plasmid design and the specific implementation of that design that we built and measured?

This sounds like a job for the CombinatorialDerivation class, which is designed for representing libraries and the relationship between derived designs and the originals.
About CombinatorialDerivation
In SBOL 3, a CombinatorialDerivation takes a Component in its template property to represent the system that you intend to vary:

Some of the Features of the template Component are designated as the variables of the library via VariableFeature child objects of the CombinatorialDerivation, which point at the variables with their variable property.
The values that can be filled in for each variable are then designated via the variant* properties of the VariableFeature.
Finally, the distribution of values is determined by the strategy property of the CombinatorialDerivation and the cardinality properties of each variable.

Representing Random Mutagenesis
In the case of random mutagenesis, if I'm understanding your design correctly, you would represent it with a CombinatorialDerivation as follows:

The template is Component of type DNA with the original sequence. Each site designated for mutagenesis can be designated with a SequenceFeature.
For each of the mutagenesis sites, create a VariableFeature in the CombinatorialDerivation with the corresponding SequenceFeature as its variable and a set of variant values for each of the mutant options set by your primers. The cardinality is probably one (if every site has to be hit by something in your protocol) or zeroOrOne (if sites might or might not be hit).
The strategy for the CombinatorialDerivation is sample, since only a few of the many possibilities will actually be realized.

Representing Samples and Characterized Constructs
Now, when you take it to the lab, you can represent your aliquot with an Implementation whose wasDerivedFrom points to the CombinatorialDerivation, representing the fact that you've got a randomly generated library whose contents are not precisely known.
Finally, once you've got your isolated products represented as Components, there is a recommended best practice to use wasDerivedFrom links to connect the derived Component and its varied Features back to the CombinatorialDerivation and its template. If you make a Collection of your results, the Collection should also link to the CombinatorialDerivation.
This should then capture the relationship between all of the products, their associated data, and the original mutagenesis design.


Generate an approapite output from DNaSP6 to Arlequin program I am performing a haplotype and nucleotide diversity of my Sanger sequences (from Watermelon mosaic virus coat protein region) by Arlequin program with the aim to see the probability that two randomly sampled alleles are different, and the expansion and genetic variation of the virus respectively.

Some of my sequences have gaps (-) because the PCR reaction could not cover all my interest region. After make the alignment in Mega and export in a Nexus file to DnaSp6 I saw that there are many option to export to Arlequin programn since DnaSP6. I do not what is the best option to generate the output from DNAsp6 to Arlequin because I must choose one option of the sites with alignment gaps option and another option from invariable sites option.

These are the options:

Sites with Alignment Gaps option:

Not considered: These sites are ignored (complete deletion).
Considered: Gaps are considered just like another nucleotide variant (fifth state).
Only gaps are considered: Only gaps information is considered to built haplotypes.

Invariable Sites option:

Removed: Invariable (monomorphic) sites will not be included in the output file.
Included: Invariable (monomorphic) sites will be included in the output file.

From my point of view I think that It is more certain for my analysis (above explained) the option of Not consider gaps (because they can buffer the information of the other sequences) and included invariable sites option because I want to analyse these paramether of my dataset.

Is it sense? I would like to know what is the best option for my dataset.

Thanks in advance.

Yes, your question makes alot of sense and its a good question because hits at the heart of differences between population genetic analysis and phylogenetic analysis. However, when you speak about gaps you are referring to 'indels' (insertion/deletion).

Indels The conservative approach is to simply delete sites with gaps. Arlequin uses a lot of maximum likelihood solutions and gaps would automatically be deleted in an Arlequin max. likelihood analysis, because it is the output per alignment per matrix. If the matrix doesn't have a 'gap'.

Indels and pop gen Thats tricky, it is better to remove them. HOWEVER the indel may contain underlying population signal, if its a maximum likelihood analysis (pop gen does this too) I suspect they will be ignored. If its pairwise, such Fst including indels might give you a better answer because there is more polymorphisms to assess population structure. So if you're doing a pairwise analysis, I'd do it both ways and then report the methodology correctly if there are differences in results and simply explain the theory.

Invariable sites This is different from an indel. An invariable site is a site which carries no phylogenetic signal, it carries information but a phylogenetic analysis could draw any conclusions because all the nucleotides/amino acids are identical at a give alignment position.

The output you obtain in allele frequency will be different and here's the rub ...

Phylogenetics analysis Arlequin is really pop gen, but it will contain some phylogenetics programs. If you see anything with 'Jukes Cantor correction', 'Kimura 2 (or 3) parameter model', 'HGW model', '"REV" model' (eek I can't remember its name any its a reverible 6-parameter model of nucleotides),  you are dealing with phylogenetics. Keeping invariable sites is normal here because the rate matrix is calculation "per alignment site".

Population genetics analysis Fst, Fis, linkage disequilibrium statistics, REMOVE THEM. Why? Because you want SNPs (single nucleotide polymorphisms). The invariable sites will get in the way of the statistics. For example is the locus recombines and you want to assess a linkage disequilibrium statistic (recombination), invariable sites will recombine, but the algorithm can't see this because the nucleotides are all the same. The only thing the algorithm can see and analyse are SNPs, which 'representative' of all the surrounding invariant sites. Thus removing invariants is essential. Arelquin might do this automatically (but I don't know the preparatory algorithms), but if you submit multiple loci Arlquin will want to know which are different loci. It will then recode the loci automatically if a single locus comprises more than one SNP. We are going into the underlying theory of linkage disequilibrium here and I'm not sure if say Fst/Fis is particularly bothered with this or is just concerned with the populations of genetic diversity.

Summary So gaps ... if they are genuine indels and not just 'N', ie indeterminant nucleotide calls, remove them for phylogenetics. If its pop gen do the calculation both ways and assess the output, if there is no difference report the 'gap/indel removal' result.

Invariant sites .. for pop gen you remove them but Arlequin might do this automatically to prepare a data set for analysis (I don't know). For phylogenetics you never remove invariant sites.

Goodluck


How can I identify a recessive and dominant gene? For example, I have two allelic genes. How can I identify a recessive and dominant gene? Are there any databases with this information? Or the answer to my question is to study the concept of "haplosufficiency"?

I'd probably look up the genes in OMIM, which provides detailed description of the functions and phenotypes associated with genes, including inheritance mode.

For example, this is what you would get for TP53:


Peak-calling using homer I have a total of 78 ATAC seq samples from which Im trying to do a peak call .I tried this batchParallel.pl findPeaks peaks function but i couldn't find the output what it makes.

So I tried using loop but again the peak file is not written due to incorrect argument usage I guess.

These are my sample directories where i have all the tags

SRR2920564  SRR2920565

I used homer to make tag directories and now i using it for peakcalling.

here is my command.

for d in */ ; do findPeaks $d -style factor -o "${d}.txt" ; done

What i need is all the peakfile should be written as the directory name.

Such as if directory is "SRR2920564" then my peak file should be written in the same directory as SRR2920564.txt

I do run the loop it runs successfully but no output written. May be a mundane mistake am not able to figure it out.

Any suggestion or help would be highly appreciated .

You just need to strip the / from $d, for which there are a number of options:

for d in */ ; do
    dname=`basename $d`
    findPeaks $d -style factor -o ${dname}.txt
done

or

for d in */ ; do
    dname=${d%%/}
    findPeaks $d -style factor -o ${dname}.txt
done

There are likely other ways one could go about this. Having said that, since you're not running these in parallel you could simply move the default output file (peak.txt) to a file name of your choice (using one of the strategies above).


Why can using more reads lead to a lower quality assembly? I am experimenting with adding additional reads to the input files I'm giving SOAPdenovo2, and there comes a point where a good contig I've been watching actually stops showing up. Does anyone have a quick answer to why that might happen?

The more reads you add the more errors you add into the assembly. This is because additional duplicate reads don't add nodes/edges to the de Bruijn graph, but those with errors do. By preselecting those reads aligning to your desired contigs, you may be further exacerbating this effect and further hindering SOAPdenovo2's ability to do k-mer correction (after all, the observed k-mer frequency of errors would be higher). This then results in additional incorrect paths through the graph that the assembler has to traverse.