Functional Annotation vs Functional Enrichment--which one should I use for Network Analysis? I am currently working on a network pharmacology project where I am supposed to find out the possible diseases our drug might have an effect on. I am not a bioinformatics student, so my knowledge in the subject matter is very limited. I have all the targets and I am supposed to get the pathways associated with the targets. I used STRING for the purpose.

I got both KEGG functional annotation and enrichment for my list of targets.
I have to build a two-mode network

Question I am not very clear about

the difference between enrichment and functional annotation and
which data should be used for the next step.

The answer to your main doubt about "annotation vs. enrichment" is that you need both, but you already have one annotation.
You can only perform a functional enrichment analysis after you have a functional annotation. In your question, you state that you have KEGG annotation for your targets. That is it, knowing which target "belongs" to each KEGG module is an example of functional annotation of your targets. Several tools can annotate genes, proteins, etc. Knowing which protein is related to each drug or disease is also a functional annotation.
Now, knowing which targets are significantly important can be achieved by functional enrichment analysis.
For the network, both information can be used together. For example, you could link all elements to their functional annotation and use the available enrichment score as a quantitative parameter of importance for calculating correlations.

I would construct two networks using both data sets and compare their congruence.

If there are major statistically differences between the two data sets that is the point at which you would ask you question.
You then assess the biological rationale for the difference in comparison to observable differences between the networks and see whether this makes biological sense.

This is a standard statistical approach, particularly when you are uncertain of the quality of the underlying analysis in relation to your biological question.


Is there any value in scaffolding the output contigs of MEGAHIT assembler given a metagenomic dataset? As far as I understood, for most assembly programs, the scaffolding step takes into consideration paired-end information in order to get from contigs (contiguous sequences) to scaffolds (longer sequences that might have some N-filled gaps).

My assembly software of choice, MEGAHIT, uses paired-end information to build the contigs, but it does not output a standard scaffold. So I am wondering the following, related things:

-is it meaningful to run a scaffold program on the output of MEGAHIT? I imagine there might be some instances in which paired-end information could span a gap.

-Which software would you recommend for it? (I've tried soapdenovo2 and SSPACE, but they appear not to be actively maintained so I have the issue of 'it doesn't work and I can't do anything about it')

-Could relevant information (eg. two contigs being connected by pair-end information) be recovered by alternative and perhaps more user-friendly means, such as exploring the assembly graph with Bandage?

Thank you for your time!

Update 2: It looks like your approach has actually been suggested here as one way to use the PE information. I guess megahit may not be really using the PE information anyways. I still believe that it's kind of weird but other people do suggest it, so maybe it's worth trying what Torsten suggests.

I think that scaffolding with PE information from the reads used as input to Megahit is somewhat sketchy. You could still try it, but I'd be worried about artifacts, just because de novo assemblers are so heuristic.

However, I think that it is perfectly ok to scaffold using orthogonal data. Here are some examples of orthogonal data:

A reference assembly (weird for metagenomics but still possibly ok).
long reads, e.g. ONT or PacBio (low coverage probably fine).
Proximity ligation data, e.g. Hi-C reads.
Mate-pair reads (also lead to artifacts at a somewhat high rate, but technically valid).
Linked-reads (e.g. 10X).
Probably other methods as well.

The tools used in each case would be somewhat different. For a little more information about how you might use this, here is a recent review of tech.

Full disclosure: I work for a company that sells Hi-C kits for such applications.

Update: Realized that I missed one part of the Q. I think that visually exploring the assembly using e.g. bandage is always a good idea. Quite possibly you can make some scaffolding decisions that way, but to me it sounds somewhat painful to do, especially in a metagenome where there are going to be a lot of multiple-branching collapsed regions.

You can use SOAPdenovo-Fusion to scaffold contigs produced by MEGAHIT as suggested by one of the developers: https://github.com/aquaskyline/SOAPdenovo2


How to import data from cell ranger to R (Seurat)? I will have some scRNA-seq data. The goal of the experiment will be to see if there is any difference in gene expression between treatment groups using the package Seurat from R.
I have read a tutorial how to do the analyze, but this tutorial does not explain how to import data.
I am a total beginner with R, I searched the internet,but the majority of the explanations related to csv/excel/txt... common files.
Data are from Cell ranger and spread in 3 files with following file extensions : .tsv and .mtx (barcodes.tsv, genes.tsv and matrix.mtx).
Collaborators ran Cell Ranger and gave these cell ranger output files : barcodes.tsv, genes.tsv and matrix.mtx.
Can someone give me the code to import these kind of data to R ?
Thank you in advance.
update: By "directory" you means "folder" on my computer for example ? So I need to save the 3 files under one folder and then to put name of the folder into brackets after the Read10x function ? I am a total beginner with R (I have used to R commander for statistical tests but that's it...).
update 2: I successed to import data into Seurat with the explanation of StupidWolf. I do not know C language ( I do not know almost any language ..). My system crashed a little further in the analyze after importing, not enough RAM... so now I am learning Linux commands to switch on the cluster.

Have a look here and the entire tutorial:
https://github.com/hbctraining/scRNA-seq/blob/master/lessons/03_SC_quality_control-setup.md

To convert data from one format into another, you may find some software tool sometimes. But mostly no such tool can be found to satisfy all the requirements that we want.
In such situations, I recommend that you write a C program that reads data from the input file and writes it into the output file. In this way, you have complete flexibility so that you can determine the exact output format and what part of data to ignore or include. To do that, you should be familiar with some concepts in C language such as arrays, loops, strings, and file operations.
For an example of such conversions, you can get ideas from the following page:
https://stackoverflow.com/questions/55347322/c-program-to-convert-a-text-file-into-a-csv-file
Of course, such a converter program can be written in other languages too such as C# and Python.
The program can be developed in R language too. You need to be familiar with import functions, loops, strings, and arrays in this language.
We cannot expect that custom tools exist for all our file convertions. Every bioinformatician must be skilled in a standard programming language such as C, C#, Java, Python, or R. Then, he/she will be able to develop any custom program that reads an input file and generates a custom output file.


How to find possible viral subsequences in a Eukaryote genome? I have an assembly of an algae and want to find if it has anything analogous to endogenous retroviruses' subsequences in its DNA.
My guess is I should manually blastn all subsequences until I find (or not something) that could have come from a virus.
Is there a better way to do this?

Your guess sounds reasonable to me but you should consider using KRAKEN2. You can create your custom database (you would have done it with makeblastdb anyway) and kraken will check k-mers against the database to find homologous sequences.

Searching directly for a reverse transcriptase may find recently incorporated retroviruses.


What tools can I use to look up my alleles, genotypes, or phenotypes from my sequenced DNA (WGS)? (I am not a Bioinformatics expert, please forgive and educate me if I've used any wrong terms or assumptions here)
I bought a "Whole Genome Sequence" kit, which gave me the following VCF files (this service did not include any analysis, just the raw data):
cnv.vcf
indel.vcf
snp.vcf
sv.vcf
My understanding is that these files show the differences of my DNA relative to a reference human genome (hg19 in my case)
What tools and resources can I use to find out which known alleles/variants of genes I have?
An example: Blood type
According to Wikipedia, there are 6 common alleles of the ABO gene, and which alleles I have determine my blood type.
Looking at my VCF results, I have these SNPs and INDELs within that gene:
<removed>

Given this, how would I figure out my blood type?
I understand that my data may not be 100% accurate, and that many phenotypes are the result of more than one gene, environmental factors, or are not well studied at all, but this should work for some simple things ie blood type, COMT, BHLHE41, CYP1A2, etc

If this is your own DNA and you are trying to analyse yourself you should talk to a genetic counselor. I personally think it's cool, but I deeply understand the potential consequences.
What you want is a variant annotation tool. There are many, but I normally write my own. Which is simply a case of pulling down a database and using pandas to match variants with annotation.
There are many databases. The biggest popular database is gnomad and the biggest disease database is clinvar. Others include omim, snpdb and 1k genomes. I'm not sure about blood sorry but I'm sure you can find one.


Fast and reliable alternatives to blast After some unexpected results (and previously reported) I heard that there are other tools for finding similar sequences besides blast that are faster and more accurate.

I only found hmmer, but I don't know if it is faster and more accurate.

Is there any review comparing such tools? A quick google scholar search didn't report anything useful.

I want to search with short DNA sequences (~400), mainly interested for taxonomic imputation so I'll be looking for similar sequences and not distant homologs. I won't try to model genes and the sequence will be unspliced as it is from direct sequencing of a a region of the 16S gene.

Besides generic nucleotide aligners, there are also more specialized tools for the alignment of 16S amplicon sequences, e.g. SINA (article, software) which is part of SILVA

You can give a try to this tool: dnaservic.es
It's actually available experimentaly. It is a real local aligner, not just a search tool like BLAST. It really returns an alignment. It is usually slower than BLAST, but bot much slower in most cases. In contrast it usually returns alignments with better scores and not only partial alignments.


Error after trimming illumina adapters I am removing illumina adapters of the NGS data with a loop. My NGS data is storage in /data/HTS_seq/.
I used this function:
 for infile in /data/HTS_seq/*_read1.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

 for infile in /data/HTS_seq/*_read2.fastq; do outfile=$infile\_trim.fastq cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile; done

and It gives to me this error:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'
  Traceback (most recent call last):
  File "/usr/local/bin/cutadapt", line 10, in <module>
cutadapt.main()
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 652, in main
options.discard_untrimmed)
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/scripts/cutadapt.py", line 327, in trimmed_and_untrimmed_files
trimmed = untrimmed = xopen(output_path, 'w')
  File "/usr/local/lib/python2.7/dist-packages/cutadapt/xopen.py", line 178, in xopen
return open(filename, mode)
      IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which is my error?
Thank you in advance.

First of all, you aren't actually setting the $outfile variable. Your loop has this basic format:
for var in *; do newVar="something" command $newVar; done

For example:
$ for var in 1 2 3 4; do newVar="new:$var" echo "$var: $newVar"; done
1:
2:
3:
4:

As you can see, newVar is never actually being set. That's because setting a variable is its own operation. So you need to do that, end it with a control character like ;, or &, or a newline, and then run your cutadapt:
$ for var in 1 2 3 4; do newVar="new:$var"; echo "$var: $newVar"; done
1: new:1
2: new:2
3: new:3
4: new:4

Or, with your data:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="$infile"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

But that would give a weird file name, like foo_S1_L001_read1.fastq_trim.fastq. You probably want names like foo_S1_L001_read1_trim.fastq instead, so you would want:
for infile in /data/HTS_seq/*_read1.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
      -e 0.1 -O 3 -m 30 --max-n 2 -o "$outfile" "$infile"
done

for infile in /data/HTS_seq/*_read2.fastq; do
    outfile="${infile%.fastq}"_trim.fastq
    cutadapt -a AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTA \
      -e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Now, the specific error you show makes me think that isn't the command you ran. Could you perhaps have run this one instead by mistake (note the / instead of \ before the _trim.fastq)?
for infile in /data/HTS_seq/*_read1.fastq; do
 outfile=$infile/_trim.fastq cutadapt -a AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC \
-e 0.1 -O 3 -m 30 --max-n 2 -o $outfile $infile
done

Your error message is:
  IOError: [Errno 20] Not a directory: '/data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq'

Which means that cutadapt tried to open /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/ as a directory and then look for a file named _trim.fastq inside it. The only way I can think that this would happen is if you ran the command I gave above and not the one you show in your question.

There is a bug here
/data/HTS_seq/*_read1.fastq

This will go to your root directory / and there is definiately no folder called data in the root. To prove it type ls / and that is where the command is pointing, then compare this to ls ./.
To solve it go to the directory HTS_seq and type pwd then copy and paste the directory structure. I suspect ./data/HTS_seq will solve it
 for infile in ./data/HTS_seq; ...

This will look for the folder data in your current working directory rather than the root directory.

Okay you've been assigned /data as root ... no biggie.
PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq

Why is / in the file name? What happens if you ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq/_trim.fastq?
and ls /data/HTS_seq/ Finally could you please ls /data/HTS_seq/PV589_03664AAB_AGGACT_read1.fastq\/_trim.fastq?
The \/ is important. What I think has happened is a / is in the file name and that will cause problems because it will think that _trim.fastq is a directory when its part of a file name, and would explain the error.
Then do the same thing using head instead of ls see what happens.
If the last ls and head works then rename the files by typing rename in the command line. If rename works it will be easy to remove / from all your files using the power of the rename function (its a Perl function BTW).


Chromosome information for miRNA Apologies if this question is previously asked. I have a BAM file alignment to hg19 and miRNA gff3 file from mirbase for annotation of BAM file. I want to extract information in the following format (e.g.):
miRNA        chr    chr_start    chr_end    expression_count    sequence
hsa-mir-155  chr2       ...         ....         1020            ATGC....

I know that I can use htseq-count for BAM annotation and for expression counts but how can I get their chromosomal location. The reason for doing this specific way is that I want to show them in a CIRCOS plot.
Just counting is not the goal, their chromosomal location also needs to be extracted. So just subread/htseq-count is not sufficient for extracting genomic location information.
Previously, I had the same idea about processing the sam records. If I want to annotate chromosomal information with miRNA using .gff file, I have to check for each line in .sam records whether their chr_start and chr_end lie in the bandwidth of start and end of any miRNA in .gff file. This process can be very heavy and slow because .sam records are of more than 5 GB. That's why I want to know if there is any extremely fast and generalized way to do it.
I want to avoid processing sam records in python because my experience with extremely heavy data like fastq or sam in python is not good. I would rather prefer awk/bash or any low-level language.
Thanks.

you can do samtools view pathtobamfile, this will convert the file to .sam format and send it to stdout, the columns can then be parsed out which contain the information you are after.

The gff should have the coordinates for each feature.  I'd use featureCounts from subread to count up how many reads hit each feature.


Automating download of results of Drug and Drug Target Mapping The drug and drug target mapping tool shows the top 3 PDB IDs for each drug based on the drug target sequence similarity search and ligand id search. Is there a way to get all the PDB ids of the drug-target PDB files and not only the top 3?

I know you can get this information after clicking the "find" button under the last column of the output table (the column named "All PDB"), but I have a query of several hundreds of drugs and I'm not willing to do it manually for each result.

Simpler option would be to not use the query interface  and directly parse the CSV file found on that page.


Extracting modules from correlation and then module clustering This is the paper where they have first done co-expression analysis of TF and then found coexpressed modules which were further subjected to module clustering ,paper

So now im clear about the first part of doing the coexpression ,the part which they have written on the paper Im quoting

  Using Pearson correlation coefficient of 0.4 as a cutoff to define a
  positive correlation between TF pairs, we obtained 37 tightly
  related TF modules that showed excellent co-expression patterns (Fig.
  4a and Supplementary Data 4).We summed the average intensity of each
  module and then performed hierarchical clustering to reveal the
  relationship among the 37 TF modules

As of now If i have done something like this ,i use all the modules coming up with no threshold as such which they have written 0.4, if i get it, then the first figure is a gene-gene pair correlation, how to filter it ,because if i get try to set a threshold i will have to look for every pair irrespective of the module found ,

How to do it in R im not sure about it, any suggestion would be really helpful,next one is how are they going for the next one which is module clustering ,as there are 37 TF modules ,each of them having different gene distribution ,here the second figure the modules are in the rows if i understand ,so how have they done it i would like to get an idea about it.

Earlier i have asked a question , i have the data as if that can be used to show me how it can be done that would be really helpful.

Co-expression in R is doable via gglot2 at the very least... I see if I can find the code.

Essentially you are performing a co-variance analysis, however Pearson' at 0.4 is a LOW cut-off and it should be 0.7.

I've got a lot on right now, so it will take a week or so before I foward the relevant code. R will have other methods outside ggplot2 however.

I was thinking ggcorr() and ggpairs() which both in the GGally package of R-cran. The output is very cute, well its ggplot2 which is trendy.

library (ggplot2)
library (GGally)
library (dplyr)

df <- select(mydata, 1:5)
cor(df) # raw output
ggcorr(df)

output: heatmap similar to Fig. a

ggpairs(df)

output: very fancy matrix plot Its pairwise scatterplots, diagrams and numbers, which is not represented in Fig. a to c, but its a very fancy version of Fig. d ... Its simplistic analysis, but its a very cute plot.

ggcorr uses pairwise (obviously) Pearson by default, which is the most versatile.

I can definitely say that 0.7 for Pearson for a straight correlation (not a covariance analysis) is conservative low and as a reviewer I would immediately complain about anything lower. A heatmap is different depending on its interpretation.

The multivariate analysis of Pearson pairwise correlation, which is what you are requesting is a different question. UPGMA style analysis directly on the heatmap, is fine if there is "one heatspot", but otherwise has limitations IMO.


Reconstruct a protein's 3D structure from a pair-wise distance matrix What is a good method to reconstruct a protein's 3D structure from a pair-wise distance matrix of residues?

I think the question is self-explanatory. I know of this method
but it seems that it uses contact maps, which actually use a certain cut-off. I don't have a cutoff based contact matrix but a matrix with actual distances.

Finding coordinates of N points that preserve (as well as possible) pair-wise distances from an NxN matrix is called multidimensional scaling.

But I don't enough about the context of your problem to tell if/how it is applicable.


How to obtain the standard gene name I have a list of genes and proteins extracted from literature. How can I convert them into standard gene names.

Gene name

Using an Ensembl ID, Uniprot ID or NCBI ID it is much easier to look up your gene. So I would strongly reccomend using one of those as a key. Personally, I'd use Uniprot, but Ensembl is probably better if your genes are human.
One reason for this is also that different resources will use different names, especially for humans.

There are a few considerations, so I will first say how this could be done properly...

Custom Uniprot matcher with weights

Some names are used by different protein, so which one is which would depend on context, co-occurrence with another gene.

Your dataset might be extremely noisy. Therefore, you likely have assigned a confidence score to your matches, yes? It would be nice to give a score to the name based on whether it is the first name of the protein or the fourth synonym in the list.

Therefore to do a conversion properly you really ought to make your own conversion system based upon context and Uniprot data. Uniprot contains the most data on synonyms and notes and binding partners etc. The dataset/API is in XML and Biopython does not have a parser, but I have some notes on the matter and can share some code if needed as I have done something similar (e.g. JSON for human gene/protein names to Uniprot IDs).

Lazy match

The NCBI entrez API is rather flexible and good at presenting to you matches for a strange names. So I would suggest using that.

bioDBnet might have what you're looking for if the gene name is in standard format

https://biodbnet-abcc.ncifcrf.gov/db/db2db.php

For proteins Uniprot can be used to look up a protein and find it's information

https://www.uniprot.org/


How to import large .bed, .gff, .vcf, .paf, .sam files into an SQL database? Are there best practices to load different bioinformatics file formats such as VCF, BED, GFF, and SAM to SQL databases? I am wondering how people out there do that efficiently.
All of these three formats are tab-separated files, so basically the following should work. I feel weird about it since most people I know don't use MySQL to work with these files.
LOAD DATA LOCAL INFILE 'bed.bed' INTO TABLE bed-file FIELDS TERMINATED BY '\t' LINES TERMINATED BY '\n' IGNORE 1 ROWS (list of the columns) SET creation_date  = STR_TO_DATE(@creation_date, '%m/%d/%y');

Answer from @liam-mcintyre converted from comment:
I don't use dask as it doesn't support enough pandas functionality (unfortunately). With pandas I do it with read_csv... if its big then read in chunks and send chunks to separate threads. If you want to ask a specific question with example data etc then I can show code.


Purpose of ### (3 consecutive pound signs / hashtags / octothorps) in GFF3 I downloaded the annotation of the C. elegans genome in GFF3 format from Ensembl.
I typed the following command, hoping to get the header of the file (lines starting with #).

grep '^#' Caenorhabditis_elegans.WBcel235.95.gff3

This returned 46798 lines, many many more than the dozen or so header lines. After searching the file for the # character, I noticed many lines that contain only three hashtag (#) symbols.

Why are these lines there?

The three # are used for splitting group of features that belong together, e.g. a transcript and it's exons. Sometimes you see with a blank line instead of the #.

To grep just the header, which has 1 or 2 # from the line start on, you can use extended regular expression:

$ grep -E "^#{1,2}[^#]" Caenorhabditis_elegans.WBcel235.95.gff3

Which means: "grep all lines that started with 1 or 2 # and are not followed by #"

If you like to create a new gff3 to get rid of this separator lines you can do this:

$ grep -v "^###" Caenorhabditis_elegans.WBcel235.95.gff3

The -v prints all line that doesn't match the pattern. (I guess this is what @Michael mean by "negator")

Hmmm, I dunno. If I was doing this I'd want an easy way to parse the file and ### would be a good marker. Anyway just pop a space character after the #. Using sed ...

sed -n '/^# /p' < ./Caenorhabditis_elegans.WBcel235.95.gff3

That should remove the ### lines

Also head -20 and sed -n '/1,20/p' ./Caenorhabditis_elegans.WBcel235.95.gff3 would get around the problem with a bit of tweaking.

The grep verion is (removing unnecessary jargon)

grep '^# ' Caenorhabditis_elegans.WBcel235.95.gff3

This is the same command,

grep '^#\s' Caenorhabditis_elegans.WBcel235.95.gff3

The point behind the code is there is no negation operator, in sed (I assume grep is the same). In Perl you would simple [^#]{3} (need checking) in your regex (remove all ###). In sed (or grep) you can exclude them using the ^ (beginning of a line) and suffixed with a "space" with only a single character in the middle all 3 letter patterns are removed.

Another approach without the space character is

sed -n '/^###/d; /^#/p' ./Caenorhabditis_elegans.WBcel235.95.gff3

Albeit this may require,

sed -i '/^###/d; /^#/p' ./file


Identifying outlier loci in R I have done some genome scan analysis and I want to identify the outlier loci. For an expected value of the FRD equal to q = 10, a list of candidate loci can be obtained by using the Benjamin-Hochberg procedure as follows:

q=0.1
L = length(adjusted.p.values)
W = which(sort(adjusted.p.values) < q * (1:L) / L)
candidates = order (adjusted.p.values)[w]

and this is part of the outliers and my output

head(candidates)
 [1]  3  4  9   10  13

But I need the output in a different type.
I want to assign TRUE if a locus has been chosen as outlier and FALSE if it has not

  head (desired_output)

  [1] FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE TRUE

Any idea how can I modify the Code?

As Devon said, removing the which will give you W as a logical vector that can be used to index the sorted p-values. To get a logical vector that can be used to subset the original (unsorted) values you’ll need to go a step further:

candidates = is.element(1 : L, order(adjusted.p.values)[W])

However, I’m not sure this will work with your data: the name adjusted.p.values suggests that these values are already FDR-adjusted. So you cannot/should not perform additional correction on them. If they are already adjusted, then you can get candidate loci simply by

candidates = adjusted.p.values < q

Remove the which(), which is converting the boolean vector to a vector if indices to TRUE.


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

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

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

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

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

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


tool to visualize a collection of phylogenetic trees I have a collection of phylogenetic trees produced from clustering DNA sequences, mostly a few nodes in each cluster/tree, with several clusters of size 1. And I would like to visualize them all together to get an idea of the relative sizes of the trees in the collection of clusters (cluster set).

E.g.

(ABC:0.0,('CL-118070':0.02010925,'CL-117615':0.00526463):0.00443474);
(DEF:0.0,('CL-117360':0.00967856,'CL-117181':0.00967856):0.00990975);
(GHI:0.0,('CL-123443':0.00123443,'CL-123443':0.00132446):0.00123435);
(JKL:0.0,('CL-134123':0.00123441,'CL-545443':0.00343446):0.00134335);
((((CL-114948:0.00908659,(CL-114989:0.00324875,CL-114959:0.00536289):0.00037038):0.00457230,(CL-114851:0.01291841,CL-114831:0.00590321):0.00077034):0.03411459,XYZ:0.00251364):0.00207594,(CL-114950:0.03074481,CL-114699:0.02256162):0.00733986,CL-116426:0.01643848);

I am thinking there must be a tool where I give it an input file with a tree per line and shows them as all in one go.

Any suggestions:

EDIT: looking at the options below:

https://en.wikipedia.org/wiki/List_of_phylogenetic_tree_visualization_software

FigTree: in this tool, one can load several trees, one per line, and there is a 'Prev/Next' menu to iterate over them, but no apparent option to display all of them at once.

Phylo.io: in this tool, one can load several trees, and display one at a time or two of them, with colouring options that compare the common nodes/branches. There isn't an apparent option to visualize many trees in one go.

PhyBin: this tool (https://cs.indiana.edu/~rrnewton/projects/phybin/) seems to cluster the sets of trees, but not visualize all of them together.

GrapeTree: this tool () can visualize only one tree at a time. It displays the last tree in the list if more than one are copy+pasted in the input text box.

There is another group of tools under Graph drawing:
https://en.wikipedia.org/wiki/Graph_drawing

Cytoscape: this tool is designed for network visualizations, which could be an option. Unsure if it visualizes newick trees, but there is are tools to convert trees to networks that could potentially work in combination with Cytoscape, e.g. https://github.com/MonsantoCo/Newick2Network

I recommend using densitree: https://www.cs.auckland.ac.nz/~remco/DensiTree/

And rwty: https://github.com/danlwarren/RWTY

Densitree plots the superimposed set of trees in a heatmap-like format:

And rwty, a very powerful R package that allows visualisation of treespace in an NMDS format (among other methods):

Dendroscope is worth a look.  I would also look into the R phangorn package.


Extracting specific sections from FASTQ determined by positions Given a FASTQ file, I'd like to generate a new FASTQ file including only certain subsections of the original sequences specified by their positions.

For example, say I want to extract the nucleotides at position 1, 3, 5-7 from each sequence:

Input:

@id1
AACCGGTCC
+
123456789
@id2
TACCGGCCC
+
123456789

Output:

@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

Note, the phred-scores are constructed such that they correspond to the position for easier reading.
Note2, the phred-score subsections should be included accordingly (i.e. the same positions)

Ideally, a command should work this way:

subselect 1,3,5-7 input.fq > output.fq

Having looked through several poplar software, I couldn't find anything working out of the box. I was thinking to use either awk or even cutadapt (systematically trimming from the 5` and pasting) but maybe out of an box solution does already exist which evaded me.

For the specific regions you mentioned, you could do this in awk. It won't be fast for a large file, but it will work:

$ awk 'NR%4==2 || NR%4==0{$0=substr($0,1,1)substr($0,3,1)substr($0,5,3)}1 ' file.fq
@id1
ACGGT
+
13567
@id2
TCGGC
+
13567

$0 is the current line. The script above will replace the line with three substrings of itself, one starting at position 1 with a length of one, one starting at position 3 with a length of one and the last starting at plsition 5 with a length of 3. The NR%4==2 || NR%4==0 ensures that this substitution only occurs on lines whose line number modulo 4 is either 2 or 0, so every 2nd and 4th line of the file.

IMPORTANT: this assumes you always have only one line of sequence in each section. The fastq format allows multiple lines of sequence (and quality scores). So your data might not conform to this. However, for regular short-read data, it should be fine.

And here's a slightly more sophisticated version in Perl, with some rudimentary error checking and which can take regions in the format you asked for (1,3,5-7):

#!/usr/bin/env perl

if ($#ARGV != 1) {
  die "Need exactly 2 arguments.\n";
}
my $file = $ARGV[0];

my @regions = split(/,/,$ARGV[1]);
my (@startPositions, @lengths);
for my $arg (0..$#regions) {
  if ($regions[$arg] =~ /^\d+$/) {
    push @startPositions, $regions[$arg];
  }
  elsif ($regions[$arg] =~ /^(\d+)-(\d+)$/) {
    for my $i ($1..$2) {
        push @startPositions, $i;
    }
  }
}

open(my $inputFileHandle, '<', $file) or
    die "Failed to open $file for reading: $!\n";

while (<$inputFileHandle>) {
  if ($. % 4 == 2 || $. % 4 == 0) {
    my $line = "";
    foreach my $start (@startPositions) {
      $line .= substr($_, $start-1, 1)
    }
    print "$line\n"
  }
  else {
    print
  }
  ## If we find a line that consists only of '+' but its line number
  ## modulo 4 isn't 3, that means that there was at least one multi-line
  ## sequence that we cannot handle correctly.
  if (/^\+$/ && $. % 4 != 3) {
    die "Line $. is just a '+', multi-line sequence found in file.\n"
  }
}

Save that as foo.pl and then run like this:

perl foo.pl file.fq 1,3,5-7

This looks like a job for samtools fqidx (http://www.htslib.org/doc/samtools-fqidx.html), assuming that your file is not enormous... read the warning in the link above.

usage would be

samtools fqidx yourfile.fq <region>

where region is specified using the fastq record identifier and coordinates.


Biplot scores from canonical correspondence analysis I'm using the R package vegan to perform canonical correspondence analysis (CCA). As input we have two matrices, one being (sites)x(species) and the other being (sites)x(conditions).

Sample data (and source of plot) are here.

Species loadings are easily accessed with summary(cca_model)$species. What I'm trying to find is the loadings for the explanatory variables, the conditions. The only summary I can find is biplot scores. Looking through the documentation for vegan I can't find any description of how they are computed. Can I sum them across CCA components to get an idea of how much they influence the data?

This is a biplot of two CCA components. The scores are used as coordinates for the arrows.

So two questions:

  What are biplot scores in the context of CCA?

and

  Can biplot scores be used to determine how much of an effect conditions have on the response variables?

Canonical correlation analysis (CCA) not correspondence canonical analysis. Biplots are the same as in PCA: the representation in a new space of the samples and the relative position of the variables on the new space.

However, CCA seeks the values of one block that maximize the correlation with the other block. Here you don't have an independent variable that you maximize from the variables you introduced on the CCA. So the interpretation is much trickier as there isn't a direction: WaterC and FallenTwigs show the opposite effect than BarSand and CoverMoss, but there is a group of variables that explain what is the effect of water on the plants(?)


wget for links inside html pages I am trying to download a file from the following repository: https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474

As you can see, there are several layers to the webpage. For example, clicking on the download tab doesn't change the URL and, the link is not a 'download link' per ce - where simply clicking the link automates a download. I have tried some of the answers on the forum, where they have advised using quotes and the operation:
wget "url/?target=link". This however does not work in the following instance.

Basically, I want to download the file labeled 'P1TLH.bam', in the download tab of the link provided:

If you want to download a file, you need to use the link to that file. Your original attempt, wget  https://trace.ncbi.nlm.nih.gov/Traces/sra/?run=SRR7276474 wouldn't work since that's a link to the trace page of the relevant run. If you want to download something else, just right click on the link (the one in the screenshot in your question), copy the URL and use that:

wget https://sra-download.ncbi.nlm.nih.gov/traces/sra64/SRZ/007276/SRR7276474/P1TLH.bam

NCBI has its own digger for this purpose. No, my bad the traditional Perl approach of LWP and HTTP::Request::Common to interrogate NCBI, isn't going to work no how much data munging or spider crawling is done, because its just an NGS data dump. ... well its diggable, but it would be easier to dump the data and  interrogate it locally.


Download data from the Human Microbiome Project via ascp I have asked this question in biostars, but I am trying here as well as people working with "omics" data might be able to help. I think my issue relates understanding how large data storage on online server works. I am trying to download large data from the Human Microbiome Project. These are available on online server and can be downloaded via aspera (ascp).

I can't figure out the right path to get the data. On their website, they give an example of command to retreive an example dataset:

 ascp -l 100M -QT user@aspera.hmpdacc.org:/ReferenceGenomes/test_genome.tar /path/to/local/directory

I executed this, using my own username, from this I am asked my password, and I can authenticate. But then I get the following error:

(Error: Server aborted session: No such file or directory)

Then, I downloaded the manifest file via the HMP portal. In this manifest, they provide the path to which these files are supposed to be available. I tried to execute again the command with one of those path, provided by the manifest file:

ascp -l 100M -QT user@aspera.hmpdacc.org:/files/5a950f27980b5d93e4c16da124000378  /path/to/local/directory

I again get the same error of "No such file or directory".

The full url provided by the manifest

I tried putting different chunks of this address into to ascp command but nothing works. I think there is something I don't understand about how to access the server via aspera and what is expected to go in this path. When I copy paste this url in google, this directly download the file I want. So I'm close, but I need to do it via command line as I need to download hundreds of these files and they are huge (6-10GB each).

How can I download the data of the human microbiome project using ascp?

I can't test for ascp because they do not longer support my current setup but you can follow the FTP instructions. With FileZilla:

url: public-ftp.ihmpdcc.org


KEGG FTP vs KEGG API I was reading the KEGG plea and I found that it doesn't forbid using the KEGG API. Then, what is in the FTP server license for personal use/academic use that it is not covered by the API?

Or I could download all the database via the API?

PS: I couldn't find how to download which genes are in which pathways but it might be I didn't expend enough time with the queries.

I found the way to download the genes and pathways thanks to this answer using KEGGREST.

It seems that the API allows to fully download the database because it allows to download info, list, find, get, conv, link, ddi.

In fact in the message of 2011-05-21 it mentions that the API includes more options to download the KGML files that previously where only possible to get via the FTP site.

We implemented Bio2BEL KEGG package to download the KEGG database via its RESTful API


Un-normalize DESeq2 counts I obtained a matrix of RNA-seq count data that has been normalized by DESeq2's median of ratio method.
I know that DESeq2 wants to take in un-normalized counts, but I do not have access to those data.
How do I best proceed here if I want to perform DEG analysis using DESeq2?
I know I can always start from .fastq files, but that would be so much extra work.
I don't think I can un-normalize the data because I don't know the scaling factors.
Any suggestions?

You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?
It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.
I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)

Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:
estimateSizeFactorsForMatrix(round(ncounts))
 sample1  sample2  sample3  sample4  sample5  sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006

If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.
Below you can use the normalized counts and you can see the difference is not huge:
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)

Compare the new count table:
head(counts(dds,normalize=TRUE))
      sample1 sample2 sample3 sample4 sample5 sample6
gene1       0       2       0       1       0       0
gene2       9       8      16       6       1       2
gene3      23      40       8      38      11      41
gene4       0       0       0       6       8       0
gene5       3       2       0       5       1       2
gene6       0       1       3      19       4       9

And your input:
head(ncounts)
        sample1    sample2   sample3    sample4   sample5   sample6
gene1  0.000000  2.2658464  0.000000  0.9310803  0.000000  0.000000
gene2  8.693480  8.3081033 16.251453  5.5864817  1.121781  2.086619
gene3 23.472395 40.0299524  7.584012 38.1742919 11.217808 40.689076
gene4  0.000000  0.0000000  0.000000  5.5864817  7.852465  0.000000
gene5  3.477392  1.5105642  0.000000  4.6554014  1.121781  2.086619
gene6  0.000000  0.7552821  3.250291 18.6216058  4.487123  9.389787

We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))

Again, I would say I agree with @ATpoint comments that if you can, start from the fastq again.. don't work with data you have doubts about.

You might be able to cheat by giving DESeq a custom size factor matrix where all the values are 1.  Then the software will feel like it's done the normalization step, but the 'normalized' values will be what you want.


Books on application of Machine Learning in Bioinformatics I'm looking for books on applications of machine learning and statistical data mining in bioinformatics with example codes provided in R and/or python.
I came across the following two books: the first being a quite good one, provides the code in MATLAB, which I'm not very familiar with. The second one too is a good book with examples in R.

Feature Selection and Ensemble Methods for Bioinformatics by Oleg Okun

Data Science and Predictive Analytics: Biomedical and Health Applications using R by Ivo Dinov

Can people suggest other books based on their first hand experiences?
Thanks

One could read Bioinformatics: The machine learning approach as an overview of the established applications of machine learning to sequence alignment, RNA structure prediction, etc...
...Yet, this book is 20 years old, and there have been many important developments in bioinformatics and machine learning in this period, which have not yet been covered in any books. In other words, reading books on bioinformatics and/or machine learning is good for getting firm basic background in either of them, but any such books are inevitably out of date at the time of their publication.
I would like however to stress the importance of the basic background - one regularly sees announcements about breakthroughs achieved by reinventing techniques that were already described in the book that I mentioned.

There are plenty of books when we talk about bioinformatics but only a few of them are available for application of machine learning in bioinformatics. From these selected books I would suggest that you go ahead with Statistical Modelling and Machine Learning Principles for Bioinformatics Techniques, Tools, and Applications. Its a recent book published in 2020 and covers most of the important topics in bioinformatics along with their ML applications. Apart from that you can also refer Machine Learning in Bioinformatics it’s quite old but still it covers most of the concepts with in-depth explanation. Along with that you can refer to this amazing book Computational Genomics with R available free online. This will cover both theory and code for different ML algorithms (both supervised and unsupervised) with their application in bioinformatics. It also covers important topics like RNA Seq Analysis, ChIP Seq Analysis in depth with code explanation of different real world examples which gives you head start in this beautiful subject.


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

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

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

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

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

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

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

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

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

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

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

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

Edit

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

I used readr for file IO:

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

Here are the timings:

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

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

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

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

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

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

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

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


Why does the SARS-Cov2 coronavirus genome end in aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa (33 a's)? The SARS-Cov2 coronavirus's genome was released, and is now available on Genbank.  Looking at it...

    1 attaaaggtt tataccttcc caggtaacaa accaaccaac tttcgatctc ttgtagatct
   61 gttctctaaa cgaactttaa aatctgtgtg gctgtcactc ggctgcatgc ttagtgcact
  121 cacgcagtat aattaataac taattactgt cgttgacagg acacgagtaa ctcgtctatc
  ...
29761 acagtgaaca atgctaggga gagctgccta tatggaagag ccctaatgtg taaaattaat
29821 tttagtagtg ctatccccat gtgattttaa tagcttctta ggagaatgac aaaaaaaaaa
29881 aaaaaaaaaa aaaaaaaaaa aaa

  Wuhan seafood market pneumonia virus isolate Wuhan-Hu-1, complete genome, Genbank

Geeze, that's a lot of a nucleotides---I don't think that's just random.  I would guess that it's either an artifact of the sequencing process, or there is some underlying biological reason.

Question: Why does the SARS-Cov2 coronavirus genome end in 33 a's?

Good observation! The 3' poly(A) tail is actually a very common feature of positive-strand RNA viruses, including coronaviruses and picornaviruses.
For coronaviruses in particular, we know that the poly(A) tail is required for replication, functioning in conjunction with the 3' untranslated region (UTR) as a cis-acting signal for negative strand synthesis and attachment to the ribosome during translation. Mutants lacking the poly(A) tail are severely compromised in replication. Jeannie Spagnolo and Brenda Hogue report:

The 3′ poly (A) tail plays an important, but as yet undefined role in Coronavirus genome replication. To further examine the requirement for the Coronavirus poly(A) tail, we created truncated poly(A) mutant defective interfering (DI) RNAs and observed the effects on replication. Bovine Coronavirus (BCV) and mouse hepatitis Coronavirus A59 (MHV-A59) DI RNAs with tails of 5 or 10 A residues were replicated, albeit at delayed kinetics as compared to DI RNAs with wild type tail lengths (>50 A residues). A BCV DI RNA lacking a poly(A) tail was unable to replicate; however, a MHV DI lacking a tail did replicate following multiple virus passages. Poly(A) tail extension/repair was concurrent with robust replication of the tail mutants. Binding of the host factor poly(A)- binding protein (PABP) appeared to correlate with the ability of DI RNAs to be replicated. Poly(A) tail mutants that were compromised for replication, or that were unable to replicate at all exhibited less in vitro PABP interaction. The data support the importance of the poly(A) tail in Coronavirus replication and further delineate the minimal requirements for viral genome propagation.
Spagnolo J.F., Hogue B.G. (2001) Requirement of the Poly(A) Tail in Coronavirus Genome Replication. In: Lavi E., Weiss S.R., Hingley S.T. (eds) The Nidoviruses. Advances in Experimental Medicine and Biology, vol 494. Springer, Boston, MA

Yu-Hui Peng et al. also report that the length of the poly(A) tail is regulated during infection:

Similar to eukaryotic mRNA, the positive-strand coronavirus genome of ~30 kilobases is 5’-capped and 3’-polyadenylated. It has been demonstrated that the length of the coronaviral poly(A) tail is not static but regulated during infection; however, little is known regarding the factors involved in coronaviral polyadenylation and its regulation. Here, we show that during infection, the level of coronavirus poly(A) tail lengthening depends on the initial length upon infection and that the minimum length to initiate lengthening may lie between 5 and 9 nucleotides. By mutagenesis analysis, it was found that (i) the hexamer AGUAAA and poly(A) tail are two important elements responsible for synthesis of the coronavirus poly(A) tail and may function in concert to accomplish polyadenylation and (ii) the function of the hexamer AGUAAA in coronaviral polyadenylation is position dependent. Based on these findings, we propose a process for how the coronaviral poly(A) tail is synthesized and undergoes variation. Our results provide the first genetic evidence to gain insight into coronaviral polyadenylation.
Peng Y-H, Lin C-H, Lin C-N, Lo C-Y, Tsai T-L, Wu H-Y (2016) Characterization of the Role of Hexamer AGUAAA and Poly(A) Tail in Coronavirus Polyadenylation. PLoS ONE 11(10): e0165077

This builds upon prior work by Hung-Yi Wu et al, which showed that the coronaviral 3' poly(A) tail is approximately 65 nucleotides in length in both genomic and sgmRNAs at peak viral RNA synthesis, and also observed that the precise length varied throughout infection. Most interestingly, they report:

Functional analyses of poly(A) tail length on specific viral RNA species, furthermore, revealed that translation, in vivo, of RNAs with the longer poly(A) tail was enhanced over those with the shorter poly(A). Although the mechanisms by which the tail lengths vary is unknown, experimental results together suggest that the length of the poly(A) and poly(U) tails is regulated. One potential function of regulated poly(A) tail length might be that for the coronavirus genome a longer poly(A) favors translation. The regulation of coronavirus translation by poly(A) tail length resembles that during embryonal development suggesting there may be mechanistic parallels.
Wu HY, Ke TY, Liao WY, Chang NY. Regulation of coronaviral poly(A) tail length during infection. PLoS One. 2013;8(7):e70548. Published 2013 Jul 29. doi:10.1371/journal.pone.0070548

It's also worth pointing out that poly(A) tails at the 3' end of RNA are not an unusual feature of viruses. Eukaryotic mRNA almost always contains poly(A) tails, which are added post-transcriptionally in a process known as polyadenylation. It should not therefore be surprising that positive-strand RNA viruses would have poly(A) tails as well. In eukaryotic mRNA, the central sequence motif for identifying a polyadenylation region is AAUAAA, identified way back in the 1970s, with more recent research confirming its ubiquity. Proudfoot 2011 is a nice review article on poly(A) signals in eukaryotic mRNA.

Not an expert, but some searching on eukaryotic positive-strand RNA viruses seems to show that polyadenylation is not uncommon.  For example, Steil, et al., 2010.


PCA plot shows big difference but not many differentially expressed genes are found I got a PCA plot of bulk RNA-seq experiment that looks the following way:

It was generated by the following code:

pcaData <- plotPCA(rld_sva, intgroup=c("Group"), returnData=TRUE)
percentVar <- round(100 * attr(pcaData, "percentVar"))
ggplot(pcaData, aes(PC1, PC2, color=Group)) +
  geom_point(size=3) +
  xlab(paste0("PC1: ",percentVar[1],"% variance")) +
  ylab(paste0("PC2: ",percentVar[2],"% variance")) +
  coord_fixed()

First sva correction was run to correct for batch effects and then rlog transformed values were plugged in to plotPCA function.

The first issue that catches the eye is that 100% of the variance is explained by just 2 dimensions. I am not sure what can one say about the data in this case. The second issue is that I get only around 20 differentially expressed genes by using DESeq2 analysis (log2FoldChange > 1, p_adj < 0.05). I know that we can not directly state that if there is a large difference on PCA there will be present a plenty of differentially expressed genes, but why is it not the case? Simple logic tells me that pca shows the difference between the samples in their gene expression, so I would expect seeing a plenty of differentially expressed genes.

You only have 4 samples total. I think it would be difficult to not have the PCA show big differences between the groups with so few points.

On the other hand, for differential expression, it is hard to get something to be statistically significant with only 2 replicates.

While this result isnt what you expected, it communicates important information about the differences along phenotype (i assume) and across replicates. PCA plots create x and y axis that represent aggregated features that account for the greatest variation, the x axis is the direction of maximal variation represents a group of genes that account for highest aggregated variance values. You can usually dissect the PCA command to find the loading (the genes used to create the axis' arbitrary unit of measurement).  The loading of PC1 (x axis)  are genes that separate the two groups along categorical groups in an unsupervised clustering technique, this is great in an analysis with substantial samples as determined by a power analysis that includes predicted effect size which need to be done in preparation of the experiment to guide the protocol. When the loading of an axis is due to a small number of genes that suggests that the variance for those genes is very large.  Again this is great if the experiment is performed with an appropriate sample. For PC2 the next greatest direction of maximal variation separates replicates and can be interpreted as indicative of the genetic diversity across a species responsible for the unique traits of individuals. These differences are far more pronounced in group 1. The suggests that the gene expression values across replicates of group 2 will yield a significantly higher correlation coefficient than those in group 1. Unfortunately this is not an appropriate sample size for any analysis really. Even the standard 3 sample comparative analysis is a statistical nightmare. PCA plots can effectively communicate magnitude and directional cohesion (or lack of cohesion) of the salient differences between groups and samples from experiments that include measurement of features in high dimensional space which is the reason they are so prominent in bioinformatics.

This plot is actually astounding as well or it would be for adequate sample numbers because PC1 indicates that between your two phenotypes or groups the are a small number of genes that account for 70% of the variation. If you repeat this experiment with greater number of samples and it yields similar results those genes that load the PC1 axis would be prime candidates for biomkarers that differentiate whatever condition you are studying and should be in vivo validated and even analyzed for potential causality.

Part of the reason that you have so few differentially expressed genes is because of your protocol as well. SVA is a surrogate variable analysis that eliminates potential lurking variables by identifying covariates or features that may yield false positives due to uncontrolled consistencies  or biases between flowcell lanes, technician performance, time to perform the experiment. If you failed to employ SVA you would have significantly more differentially expressed genes but they would be potential false positives and should be interpreted as low confidence positives. Having performed SVA before PCA suggests that there is high confidence in the role of the several genes responsible for the greatest differences across the two conditions.
Unsupervised clustering techniques require some understanding of the basic concepts of linear algebra. MITOCW has a fantastic free course by gilbert strang that would help you to more accurately interpret these types of plots as well as many others like MDS plots, heatmaps, and dendrograms. The URL of the course and its free textbook are included below.

https://ocw.mit.edu/courses/mathematics/18-06-linear-algebra-spring-2010/

http://facultymember.iaukhsh.ac.ir/images/Uploaded_files/[Strang_G.]_Linear_algebra_and_its_applications(4)[5881001].PDF

I like Rbloggers tutorials on these subjects usually but so many good ones currently exist.

https://www.r-bloggers.com/principal-component-analysis-in-r/

I hope this helps.


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

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

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

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

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

I realized that your question is specifically gene set enrichment analysis. However, I personally have never tried any R package doing this. Although topGO may seem to be the one that may work for you? What I have been doing is that for gene enrichment analysis using Gene Ontology, simply paste and run, or DAVID when it was still updating. Sorry, I am not an expert in incorporating biological significance study in the computational result.

It seems that you are also not quite sure what tools you should use for downstream analysis. I think the following publication may be helpful.

https://www.ncbi.nlm.nih.gov/pubmed/30228881

The same primary author also published a paper titled, Cluster Headache: Comparing Clustering Tools for 10X Single Cell Sequencing Data (For some reason I could not find the link to it, will edit the answer if I find it. )

Hope that helps.


Why does samtools mpileup sometimes include ref bases (other than ',' or '.')? This is my first post here. I can think of no way of giving an easily reproducible example, as per the stack ethos. SO apologies in advance and any feedback on question format appreciated.

I have called samtools mpileup (Version: 1.9 (using htslib 1.9)) on some oxford nanopore data (mapped to GRCh38).

I am seeing some instances like below (grouped with pythons Counter), here reference is T at this position:

 Counter({',': 34, '.': 33, 'c': 9, 'g': 6, 'a': 4, 't': 2, '-:ACTC': 1, '-:ACTCGG': 1, '-:5': 1, '$': 1, '*': 1})

Why is 't': 2' in there?

flags used: -A -B -Q 1 -R -x

Why is 't': 2' in there?

How else do you think the pileup should show you that two reads (in the reverse direction) have a T in that position?


How can I take cell number into account to find total RNA expression? I'm hoping to quantitatively show differences in total RNA expression for a gene in a cluster of interest between different experimental groups. My exported average RNA values for each experimental group in my macrophage cluster are shown below:
00 d Control: 10.43
00 d PLX5622: 14.48
07 d Control: 16.08
07 d PLX5622: 17.12
28 d Control: 15.90
28 d PLX6722: 19.40
However, the violin plot for Cd68 in the macrophage cluster looks like this

VlnPlot(combine.combined, features = ("Cd68"), pt.size = 1, idents = "Macrophages", group.by = "orig.ident", split.by = NULL, assay = "RNA"):

The violin plot suggests that the Control groups actually have more total CD68 RNA in them than in the PLX5622 groups, which isn't what the average RNA values show. Any code you can share showing me how to get total RNA per group would be an amazing help!
Thanks so much!

(I currently cant comment as I used up all my rep for a bounty.)
Could you explain how you calculate the average? In you plot the averages are between an expression level of 2 and 3, while you report values over 10. Given the data you show, there is no significant difference and I would guess p-values are above 0.8!
And you mention that the violin plots suggest the control groups have higher expression.. but how do you see that? For me they are identical with minimal trends that the PLX5622 groups are higher.
And to address your question regarding showing the total: No, you would want to show the mean or median, as the total is pointless, given deviating sample sizes.
EDIT: Overall, I have a feeling you uploaded the wrong picture!
EDIT2: I think I got your point now: The difference in mean is caused by an increased number of outliers in your controls at the value 0. Despite those, the means are identical and this is what can be seen through the violin plots. You should rather try to diminish the effect of these outliers.

To me, it looks like the widest parts of the plots are higher in the treated.  I don't see why you are so sure that the violin plots are hugely different from the calculated averages.
Either way, the differences are really small.  Maybe not statistically significant.
BTW, in the future, you might want to disguise your treatment and gene of interest when on a public forum.


Recommended workflow for RNA long read Kallisto-like TPM estimation? On short reads, Kallisto/Salmon is a standard workflow for measuring RNA transcript counts in TPM. However, when I tried to Google a similar workflow for long reads, it's not clear there is a replacement.

Stringtie now does long reads too. We've found flair to be quite difficult to install and to crash often. That said, stringtie isn't very good about showing you its working (e.g. I have no idea how to get it to give the assembly it uses to figure out all the isoforms). Its values for TPM/FPKM also seem a little weird, in that they're perfectly correlated with each other, and for transcripts they're perfectly correlated with coverage too.
https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1910-1
https://github.com/gpertea/stringtie
TALON is another tool that does isoform assembly. I'm not sure whether it does quantification though.
https://github.com/mortazavilab/TALON#

While the selective-alignment algorithm that is part of salmon is designed for short reads, salmon also supports quantification via alignments in a SAM / BAM file (aligned by other tools).  You have to be somewhat careful about setting the proper command line flags, but nonetheless, the quantification should work fine.  You can, for example, find a pipeline that uses minimap2 for alignment followed by salmon for quantification of long reads here.


Extract residue sequence from pdb file (in biopython, but open to recommendations) I'm new to Biopython and I'd like to extract the sequence of residues from a pdb file.
My two questions are:

What is the simplest way to do this? (Esp. when there is more than one sequence)
and
Should I be worried about "PDBConstructionWarning: WARNING: Chain B is discontinuous"?

So far, I've obtained the residue sequence via:
p = PDBParser()
structure = p.get_structure("1ppi", "1ppi.pdb")
ppb=PPBuilder()
for pp in ppb.build_peptides(structure):
     print(pp.get_sequence())

seq = pp.get_sequence().__str__()

This seems to work well for this molecule. However, is there an easier way, especially when there is more than one sequence?
For example, I've read that one can also do
res_list = Bio.PDB.Selection.unfold_entities(structure, 'R')

but res_list is not a sequence of residues in str, and I don't know how to convert the output from res_list to a str sequence.
In addition (or perhaps because I'm going through the PPBuilder), I've recently gotten a lot of warnings of the type:
/usr/local/lib/python3.8/site-packages/Bio/PDB/StructureBuilder.py:89: PDBConstructionWarning: WARNING: Chain A is discontinuous at line..
(For example, with structure = p.get_structure("5owu", "5owu.pdb"))
I've seen a discussion about how to silence such warnings, but should I be worried about this? I've also noticed that when I get these warnings, pp builder seems to give me more sequences then are there.

There are other ways to do it (many and some might be easier), but if you want to do it with BioPython PDB module you can iterate the residues like this:
 # You can use a dict to convert three letter code to one letter code
 d3to1 = {'CYS': 'C', 'ASP': 'D', 'SER': 'S', 'GLN': 'Q', 'LYS': 'K',
 'ILE': 'I', 'PRO': 'P', 'THR': 'T', 'PHE': 'F', 'ASN': 'N',
 'GLY': 'G', 'HIS': 'H', 'LEU': 'L', 'ARG': 'R', 'TRP': 'W',
 'ALA': 'A', 'VAL':'V', 'GLU': 'E', 'TYR': 'Y', 'MET': 'M'}

 # Just an example input pdb
 record = '1pa2.pdb'

 # run parser
 parser = PDBParser(QUIET=True)
 structure = parser.get_structure('struct', record)

 # iterate each model, chain, and residue
 # printing out the sequence for each chain

 for model in structure:
     for chain in model:
         seq = []
         for residue in chain:
             seq.append(d3to1[residue.resname])
         print('>some_header\n',''.join(seq))

Not really, is very common, you can turn it off like this (taken from https://lists.open-bio.org/pipermail/biopython/2014-July/015371.html):
 import warnings
 from Bio import PDBConstructionWarning
 #your code
 with warnings.catch_warnings():
 warnings.simplefilter('ignore', PDBConstructionWarning)
 #your code which might trigger the warning
 #rest of your code here

 However, because this is so common, you can just use the
 QUIET=True option on the PDBParser object:

 from Bio.PDB.PDBParser import PDBParser
 struct = PDBParser(QUIET=True).get_structure('tmp', pdbf)

Try help(PDBParser) for more, or see:
http://biopython.org/DIST/docs/api/Bio.PDB.PDBParser%27.PDBParser-class.html#__init__

You can also use the SeqIO module from Biopython https://biopython.org/docs/1.75/api/Bio.SeqIO.PdbIO.html
and do something like this:
from Bio import SeqIO

for record in SeqIO.parse(target, "pdb-atom"):
    print(record.seq)

where target is the path to your pdb file
the parser has multiple input/output functions, and pdb-atom is the one that will read out of the coordinates itself.


Interpretation of the violin plots from sc-RNA-seq I'm confused about the meaning of the black dots and the red shape in the violin plots from the seurat tutorial:

The black dots represent the values for individual cells. The "nGene" plot (the first one) shows the number of detected genes for every cell.

The red shape shows the distribution of the data. For the "nGene" plot, you can see that the average number of genes per cell is about 900 and most of the cells have roughly around 700-1100 genes.

To me, it looks like the actual data points which are used to create the violin plot distribution.
More details about the plots can help in understanding then better.


Is the genome a programming language (i.e. LISP)? Can we analyze it with computer science theory? Question moved here from the biology stack exchange.

It looks like the genome is a kind of LISP language.

Operon structure of genome:
- https://github.com/philschatz/microbiology-book/blob/master/resources/OSC_Microbio_11_07_Operon.jpg

LISP function structure:
- http://support.ircam.fr/docs/om/om6-manual/res/listprefix.png

Nested genes:
- https://player.slideplayer.com/31/9782148/data/images/img11.jpg

Nested LISP cons cells:
- https://www.tutorialspoint.com/lisp/images/treestructure.jpg

The basic similarity I am seeing is of general structure: (func var0 var1 ... varn) and nested cons cells.

The correspondence with genome structure is the operon portion, which has:

operator = function
structural genes = variables
structural genes can be nested like cons cells

Is this a global structure for the genome, then a specific prediction is we can do something simple like counting open and close markers in the nested genes and they will be equal.

So, why would such a correspondence matter, if true?  The basic idea is that if the LISP structure is consistent, then we could maybe look for other programming constructs built on the structure, such as recursion/looping, maybe templating or even object oriented programming. Or just identify simple properties like balanced opening closing markers. The mere fact it has a LISP structure would be fascinating if true.  And since LISP is essentially an abstract syntax tree (AST), which all programming languages have, perhaps this structure I see in the genome is actually an AST, which implies an even more fundamental correspondence between the genome and programming languages.

More generally, we have a wealth of knowledge reverse engineering human programming languages.  If the genome is also a programming language similar in some respects to human programming languages, such as being a LISP variant, then we could apply this same reverse engineering knowledge base to better understand the genome.

So, first of all, is my high level observation correct?  Does the genome, at least for bacteria, have a global LISP like structure?  If yes, then if the genome is a kind of programming language, i.e. LISP, can we analyze it like software to understand its functionality?  Has anyone done this before?  Simple searching on Google Scholar turns up lots of software for analyzing genomes, but nothing about analyzing the genome as software.  It would be interesting to know whether anyone has seriously examined the question before, and either made some useful discoveries, or found out there is only a very loose analogy with programming languages, and thus knowledge regarding programming languages is not of use regarding the genome.

Hopefully this question is clear enough now, if not please ask me for further clarification.  I would like to get some sort of decent answer.

UPDATE: Someone who has put thought into "DNA is code" idea.  At a more general level than whether specific genome structures match specific programming structures.

UPDATE #2: The late Scott Federhen of NCBI wrote an article comparing genomic replication with the lambda operator.

UPDATE #3: It might be the case that LISP was inspired by the genome.

I see what the metaphor has been inspired by, where I think it got turned around, and might have a better one.

I'd like to back up a bit and establish some common ground.  It isn't wrong to expect that some aspects of genetic code and software code are similar, and while biology is usually very specific, some aspects can be seen productively from the light of computer science ideas.  For example, mass-action molecular dynamics, which are often a good model for protein interactions, are Turing-complete.  Here's a fun example where people have written a programming language out of molecular dynamics:  https://arxiv.org/abs/1809.07430

It's also true that proteins have particular substructures leading to certain functions, such as using alpha-helices to react less with the outside world.  Learning about larger substructures, called protein domains, help identify the functions of families of proteins.

The problem is that the parts of the genome which appear nested like LISP cells actually aren't included in the protein.  The regulatory regions, promoters, inhibitors, introns, and so forth never become part of the key "computational" machinery.  This nesting isn't actively part of the computational activity proteins do perform.  The main quality of LISP is that the structure and function are perfectly self-similar throughout an expression, but the structure of the regulatory regions are neither replicated down into the protein nor are necessarily recapitulated on a larger scale.

Unfortunately, we don't expect the computational structure of proteins to be the same from one protein to another.  This is because "binds only with X" and "binds only with Y" must be implemented differently to successfully be distinct operations, and it is these binding constraints by which molecular dynamics implement "computational" operations, such as "not X" and "not Y", and so those are distinct as well.

I think a better way to think about the regular structures of the regulatory apparatus is more like a network protocol.  Along an otherwise undifferentiated binary channel, one has to identify when and where the data is delivered to an application.  Similarly, along a strand of bases, the machinery of translation has to have an indication of where the "application" of the protein starts and ends, and when that protein should be delivered to the cell.  We might compare some of the images of genetic structure that provided in the original question with layouts of network headers: http://slideplayer.com/slide/4798836/15/images/16/TCP/IP+Packet+Structures.jpg  We might compare bacterial versus multi-cellular regulation as having different delivery mechanisms and data structures for delivering the application information, for different functional consequences, in the same way TCP and UDP have different structures for different trade-offs.

This is all complicated somewhat by the fact that proteins do then interact with the regulation of genes, leading to interesting "computational" dynamics such as feed-forward regulation.  However, I'm not aware of any particular structure of the regulatory apparatus that necessarily correlates with these larger functional patterns.  The opposite is frequently biologically useful: the same structure of a protein helping to activate one gene might simultaneously inhibit another.

One book I've found very accessible and useful in understanding where computation is useful biologically is Dennis Bray's "Wetware: a computer in every living cell".  It gave me a clearer picture of how computation occurs in cells.  That should open up some introductions to systems biology which explains the biological-role of larger molecular "computational" patterns.

Believe it or not, there are people who know more than you do about programming, and who also know a hell of a lot more than you do about the genome.  This is not an avenue those people think is worth pursuing.

DNA is a chemical, and interacts with other chemicals based on its shape.  There is no equivalent to that in a human designed programming language.


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

I'd download the fasta, and use samtools, or maybe bedtools.


How can I find out if my gene of interest is duplicated? To identify homologous genes, I performed tBLASTn analyses against my genome of interest resulting in sequences of a closely related species. Using this method, I found two (and in some cases up to five copies) of my different genes of interest in the target genome.

How can I determine if these gene copies correspond to true duplication events?

I aligned the newly found sequences to the reference genes and I found good alignments in more than one region (20% - 50% identity), but I am unsure if this is enough to consider them as duplicates.

Is there a threshold to help me determine if they are indeed duplications? How can I distinguish between a true duplication event and a pseudogene or a paralogous gene? The target genome belongs to a polyploid organism, although the genome assembly type is listed as haploid.

Any help will be greatly appreciated.

Thank you all very much for your responses, I am now beginning to understand how to find out if my genes are duplicated.

I performed some of the analyses that you suggested and I would very much appreciate your input on the results.

First, I followed @terdon's advice and I took the best hit for one of the query (source) genes in the target genome and I used this sequence as a query against the query (source) genome, i.e., I performed a reciprocal BLAST analysis (in my case it was a reciprocal BLASTp analysis). The first hit was my original query gene, which had an E-value of 0.0. For other sequences the reciprocal BLASTp hits were more variable, ranging from 3e-42 to 2e-154, but the first hit was always the query sequence. I also looked at the protein domains and they are conserved, i.e., both the query and the target sequences have the same domains. As for the exonic structure, it is fairly similar, for e.g., the query has 12 exons while the target sequence has 10. Exons are similar in length. So, in summary, when the potential duplicates from my target genome are blasted against the original source genome, the best hit is the original query sequence.

As for @zhanf's suggestion, I did a basic phylogenetic analysis using Clustal Omega (neighbour-joining tree without distance corrections) using four sequences: the original query sequence and the three gene copies (seq_5-1, seq_5-2 and seq_5-3) that I found in the target genome. This is the result:

I also compared each of the sequences against the original query:

seq_5-1 vs query: 34% identity (query length 780 aa, seq_5-1 length 400 aa)
seq_5-2 vs query: 42% identity (query length 780 aa, seq_5-2 length 580 aa)
seq_5-3 vs query: 40% identity (query length 780 aa, seq_5-3 length 490 aa)

So in summary, these sequences would appear to be paralogues, right?

Again, thank you very much for your help, I greatly appreciate it!

Leah

You might want to have a look at What is the difference between orthologs, paralogs and homologs?. Paralogs are duplicates by definition, and whether a gene is a pseudogene or still active has no bearing on whether it comes from a duplication.

In any case, to figure out what's going in in your target genome, you should use the target's genes. So, first get the best hit for each of your query genes in the target genome. Then, to investigate its evolutionary history, use that sequence (the one from the target genome) as  query against the target genome. That will let you assess how similar the target's hits are to one another and that's how you can attempt to figure out their evolutionary history.

Note that 20% identity is very low for closely related species. While there is no golden threshold, I would expect recent duplication events (and in closely related species, a duplication event that is in one genome and not the other will be recent) to produce much more similar sequences. In any case, you can't do this based on sequence identity alone.

Are the protein domains of the gene products conserved? Is the exonic structure of the gene conserved? You should also try taking the putative duplicate of geneA from genome2 and blast it against genome1 (your original source genome). Is the best hit the geneA from genome1, the query you used to find the putative duplicate in genome2? Or is there a better hit for it against another gene?

These are all questions to investigate and take into account when reconstructing the evolutionary history. You can't just look at identity and use some threshold to decide this sort of thing.

If you blast the potential dupicate from your target genome against your original source genome, is the best hit the duplicate?

I don't have enough reputation to comment but if a sequence is duplicated then there has obviously been a duplication event. Doesn't matter how scientists call it now.It can be older or more recent but it's yet a duplication.


tool for the prediction of impact of mutations on protein localization I need to know if the mutations present in the "position of interest" can have an impact on the protein localization ability. Is there any tool to make such prediction that can make my life easier?

I know the position of nuclear localization sequence in a protein and I also know it has a mutation, but I don't know how to predict the impact of this mutation on the localization of this protein ie if this mutation hindering protein's entry to the nucleus or not? Should an analysis of protein stability be useful it this case, as an alternative?

I am not interested in localization itself but the impact of mutations on the localization.

Do you have the protein sequence after the mutation available with you?
If so then maybe you can BLAST them and find out where the best hits are localized using GO annotation (cellular components). I am aware that these won't be your exact protein sequence but they can at least give you a rough idea.


Error while Importing fastq sequences to jMHC software I am trying to analyze my .fastq files (MiSeq exon reads) with the jMHC software, but output files I have obtained have no results.

I think I might be doing something wrong during the first step: "connect to database".
Because on the "import/extracting" step, I open all my .fastq files, click "Start", and it shows this message:

  "Success! Processed: 139, Imported: 0".

So I assume it is reading the files OK, but not importing them to the "database"? Could it be the extension of my files? The software says it operates on FASTA files. Are .fastq files not compatible?

Anybody that has successfully run this program could give me a hand?

Yes, use .fasta , you can try to convert your fastq files with NextGen Workbench for example.


How to simulate NGS reads, controlling sequence coverage? I have a FASTA file with 100+ sequences like this:

>Sequence1
GTGCCTATTGCTACTAAAA ...
>Sequence2
GCAATGCAAGGAAGTGATGGCGGAAATAGCGTTA
......

I also have a text file like this:

Sequence1 40
Sequence2 30
......

I would like to simulate next-generation paired-end reads for all the sequences in my FASTA file. For Sequence1, I would like to simulate at 40x coverage. For Sequence2, I would like to simulate at 30x coverage. In other words, I want to control my sequence coverage for each sequence in my simulation.

Q: What is the simplest way to do that? Any software I should use? Bioconductor?

I am working on a Illumina sequencing simulator for metagenomics: InSilicoSeq

It is still in alpha release and very experimental, but given a multi-fasta and an abundance file, it will generate reads from your input genomes with different coverages.

From the documentation:

iss generate --genomes genomes.fasta --abundance abundance_file.txt \
    --model_file HiSeq2500 --output HiSeq_reads

Where:

# multi-fasta file
>genome_A
ATGC...
>genome_B
CCGT...
...

# abundance file (total abundance must be 1!)
genome_A    0.2
genome_B    0.4
...

I didn't design it to work with coverage but rather abundance of the genome in a metagenome, so you might have to do a tiny bit of math ;)

You can split your FASTA file sequence-by-sequence using

split -a 6 -p '^>' your_file.fa seq_

and then use any existing read simulator supporting coverage (ART, DWGsim, etc.). If you want to have all the reads mixed (not ordered by the original sequence), you can use RNFtools.

Edit 1:

As @terdon pointed out, the previous command works on OS X only. An analogical one liner for Linux (but with a slightly different naming scheme using numbers rather than letters) can be

csplit -f seq_ -n 6 your_file.fa '/^>/' {*}

To make this command work also on OS X, one needs to install coreutils (e.g., using brew) and then use gcsplit instead of csplit.

Edit 2:

Once FASTA is split by sequences, the simulation becomes straightforward and many different approaches can be used. My favorite one is using GNU Parallel. Imagine that you have your coverages in a text file called covs.txt on separate lines and in the same order as the sequences in your_file.fa, e.g.,

40
30
...

Then you can simulate reads from the original sequences using DWGsim by

ls -1 seq_* | paste covs.txt - \
    | parallel -v --colsep '\t' dwgsim -1 100 -2 100 -C {1} {2} sim_{2}

and merge the obtained FASTQ files using:

cat sim_seq_*.bwa.read1.fastq > reads.1.fq
cat sim_seq_*.bwa.read2.fastq > reads.2.fq

One possible danger of this approach is that we have assumed that the number of seq_* files is the same as the number of lines in covs.txt, which might not be true (by mistake). We should check this prior to the simulation step, e.g., by:

[[ "$(ls -1 seq_* | wc -l)" == "$(cat covs.txt | wc -l)" ]] \
    || echo "Incorrent number of lines in covs.txt"

Another caveat is that the simulated reads are not in a random order (they are grouped by source sequences).