Looping over several files in bash I have a bash script:

I am wondering how I can change this script to loop over a bunch of .vcf files and give output .txt with the name of corresponding .vcf

I tried changes done in similar script in  Changing this code in a way to work for my files but I just  messed everything up

If we save the script pasted in the main post as a sh file and we have some .vcf files in a folder, by this line we can iterate over vcf files to extract what mentioned in the script returning .txt output by the name of corresponding vcf

for file in *.vcf ; do bash indel_vcf_parasing.sh $file ; done


About optimal distance (or dissimilarity) measure in proportion data My data is proportional data, but variables are correlated with each other. So I am wondering what kinds of distance or dissimilarity we should use to take correlation into account.
If we would like to use data with proportion (e.g.,total sum scaling or compositional data), what kinds of distance or dissimilarity can we use?
Bray-Curtis is one of them, but does this make sense to use arcisin or logit transformation and use Pearson correlation?
In continuous variables, we can use Pearson and Spearman distance. In our case, what kinds of distance can we use? Pearson distance cannot be applied to data with 0-1 scale. Logit transformation -> Pearson?

I would suggest reading a little on the topic; my understanding is that distance measures can be pretty robust to proportion/compositional data. Here is a relevant textbook chapter:

Other distance measures assume that the data are proportions ranging from zero to one, inclusive.

Among these are indeed Bray-Curtis as well as Jaccard distance (and others). I think that you will want to use one of the more traditional distances rather than a Pearson/Spearman, though they are getting at the same problem. A transformation might help, I suppose, but I would prefer to choose a simple but predictable distance measure rather than a correction on a possibly unsuitable distance that might have unpredictable behavior.
It's a little hard to tell from your question, but it seems that you are wanting to in some way account for indirect correlations using a distance metric. I would not suggest trying to make your distance measure do all this work. Instead, I would use a simple but clear distance method, and then use a multidimensional analysis such as principal coordinate analysis (or similar) to understand the correlation structure of your compositional variables. This is a fairly standard workflow in metagenomics, for instance, which relies heavily on high-dimensional compositional data. I would recommend looking more into metagenomics statistical applications for inspiration if needed.
The choice does also depend to a certain extent on the specific application. In some fields certain methods are familiar and trusted even if there is not a strong empirical/theoretical rule.
You've asked a lot of questions about Pearson/Spearman correlations. Much has been written on the use of correlation analysis on compositional data, most of it cautioning against doing this kind of analysis without special treatment. This leads to difficulties in building meaningful association models for different variables. One method that I'm familiar with in metagenomics is sparse correlation, e.g. SPARCC. But there are newer methods that may be better, for example here and here.


What is the best distribution to model the FPKM values from normalized RNA-Seq data? I know that the discrete raw counts from the RNA-Seq data are usually modeled by a negative binomial or a Poisson distribution, but what I am working on are the FPKM (Fragments Per Kilobase of transcript per Million mapped reads) values which represent a transformation of the raw values. They are continuous except many actual 0s. What is the probability distribution FPKM values are likely to follow, or are assumed to follow? Even after taking logarithm twice (after replacing 0s with 1, of course), the data does not follow a normal distribution since there are a lot of 0s.

FPKM/TPM values are generally log-normally distributed. Reference : Gene expression profiling in single cells from the pancreatic islets of Langerhans reveals lognormal distribution of mRNA levels


How can I get a list of papers where a particular mutation has been mentioned? I have recently tried a few tools (MutationFinder, tmVar, PubTator etc.) that extract mutations from a given text. However, for my work I need something that will take a list of mutations as an input and return a list of pubmed id in which the mutations has been found.

Example input:
mut_list = [Y702F, L483Y]

Example output:
pub_list = [12215417, 18537569]

Can anyone give any idea?

Thanks!

You could use the pubmed API to query it directly.  Typically I use Biopython for my pubmed queries.  There are some good examples in the cookbook and on other SO posts.

This example from the cookbook link above has the guts of what you need.  You will just have to iterate over your list.  Depending on the number of queries you want to do you might have issues getting throttled by NCBI

In this example, we will query PubMed for all articles having to do with orchids (see section 2.3 for our motivation). We first check how many of such articles there are:

>>> from Bio import Entrez
>>> Entrez.email = "A.N.Other@example.com"     # Always tell NCBI who you are
>>> handle = Entrez.egquery(term="orchid")
>>> record = Entrez.read(handle)
>>> for row in record["eGQueryResult"]: ...
        if row["DbName"]=="pubmed": ...
            print(row["Count"]) 463

Now we use the Bio.Entrez.efetch function to download the PubMed IDs of these 463 articles:

>>> from Bio import Entrez
>>> handle = Entrez.esearch(db="pubmed", term="orchid", retmax=463)
>>> record = Entrez.read(handle)
>>> handle.close()
>>> idlist = record["IdList"]

This returns a Python list containing all of the PubMed IDs of articles related to orchids


plink: --update-name vs. editing the BIM I am working with data from the HGDP public files. These are massive VCFs from whole genome sequencing. To begin with, I converted these to BEDs using plink (1.9). The next problem is that these files use rsid names instead of chrpos names. I would like to change to chrpos. I am wondering whether one can simply edit the BIMs or whether I have to do use plink's --update-name call. The problem with the latter is that it doesn't support duplicates. The HGDP files are full of duplicate variant names because many of the variants don't have rsid's and are just named ".":
$> plink --bfile hgdp_chr21 --update-name updateids --make-bed --out hgdp_chr21_chrpos
PLINK v1.90b6.17 64-bit (28 Apr 2020)          www.cog-genomics.org/plink/1.9/
(C) 2005-2020 Shaun Purcell, Christopher Chang   GNU General Public License v3
Logging to hgdp_chr21_chrpos.log.
Options in effect:
  --bfile hgdp_chr21
  --make-bed
  --out hgdp_chr21_chrpos
  --update-name updateids

64259 MB RAM detected; reserving 32129 MB for main workspace.
1086522 variants loaded from .bim file.
929 people (0 males, 0 females, 929 ambiguous) loaded from .fam.
Ambiguous sex IDs written to hgdp_chr21_chrpos.nosex .
Error: Duplicate ID '.'.

Thus, one gets an error when trying --update-name and I am wondering if I can simply edit the BIM files to be chrpos instead of rsid. I can't find an answer to this anywhere. It seems that if the BED files don't contain a reference to the variant IDs, but only to the row in the BIM file, then this approach should be possible.
Here's the current BIM file:
user@desktop:/media/luks8tb1/data/genomics/HGDP$ head hgdp_hg38.bim
1   rs1331781659    0   10153   G   A
1   rs1179689498    0   10163   C   T
1   rs201694901 0   10180   C   T
1   rs199706086 0   10250   C   A
1   rs111200574 0   10257   C   A
1   rs145427775 0   10291   T   C
1   rs868413313 0   10297   T   C
1   rs112750067 0   10327   C   T
1   rs1035171912    0   10330   A   C
1   rs868804735 0   10333   T   C

As can be seen, very few lines have the "." variant name. In fact, in the first 100 lines there are 0.

I would reccomend using bcftools on the original vcf files before you convert them to plink, to fill in missing IDs using the command:
bcftools annotate --set-id +'%CHROM\_%POS\_%REF\_%FIRST_ALT' file.vcf

This means you won't have any duplicate ID names and plink should be happy about that. It will preserve the positional information in the .bim file as well.


Does MinKNOW work with Mac OSX high sierra 10.13.1? I have done a compatibility check as described here on my computer to know if I can use MinKNOW on my mac laptop. It says that my OS is incompatible and that the OSX must be Yosemite or El Capitan. Did anyone manage to make MinKNOW work on High Sierra? I find it odd that they would not support more recent version of mac.

Some replies from ONT is:

  OS X is supported and you can view the installation instructions on
  our Downloads page linked below. At the moment High Sierra is not
  officially supported as we have not yet fully tested it, but other
  users in the community have reported no issues with it so far.

  The compatibility test only performs a simplistic check to see if you
  are using one of our officially supported OS X versions. Since you are
  on High Sierra, it will show a failure. That does not mean the
  software can't be installed and used. As I said, other users have
  reported no problems, but we can't officially support it if you
  encounter an issue.

So, I am going ahead and installing MinKNOW on my High Sierra.


What is cellranger doing in comparison to other methods? I've recently started working with the 10X-Genomics platform with Illumina (MiSeq and HiSeq) for single-cell RNA-Seq. I've been recommended the "cellranger" (version 2.1.0) which I understand handles the barcoding of the platform and performs processing of the data. For example:

1) convert raw Illumina .bcl files into .fastq (demultiplex)

cellranger mkfastq

2) perform alignment, barcode counting, and generate gene-barcode-martices

cellranger count

However, it is unclear which underlying bioinformatics procedure (with which parameters) are being used here. For this pipeline, I wish to know how the trimming, alignment and counts are being performed. This is important to understand whether further normalisation and quality checking is necessary. What steps would I need to do to perform a similar analysis?
If for example, I wish to document in more detail what has been performed here in a reproducible manner or reproduce the analysis pipeline on a different platform (such as Nadia) to compare them.

Ideally, I would like to know the analogous pipeline using open-source tools but at the moment I am first seeking to understand what is going on under the hood when I am running the cellranger.

cellranger mkfastq is not necessary anymore.  It used to be that the cellranger software wanted the reads to be interleaved, and you could use cellranger to do that for you if you couldn't do it yourself.  Newer versions of cellranger will take the fastq files just like Illumina's bcl2fastq makes them.

Cellranger count aligns the reads, filters away duplicates based on UMIs, tries to figure out which cell barcodes really captured cells, and tells you for each cell barcode how many reads hit each gene.  You can go through the .bams, and see which reads were counted for what gene and what barcode, based on the flags and other tags in the .bam

If you want another pipeline to compare to, you can use the McCarroll lab's DropSeq pipeline, which is a bunch of little scripts, around the output of STAR aligner (just like 10XGenomics uses).  I have put the same data sets through both 10X and Dropseq, and they return about the same thing.


Retrieving all genes for a gene ontology term I'm looking for an easy way to retrieve all the genes in a list that are associated with a certain GO term, preferably using R/Bioconductor packages. I'm not interested in under/overrepresentation or enrichment.
For instance, I want a list of all genes known to be located in 'presynaptic endosome' (GO:009883).
I tried the method referred to in an older post (Retrieve all genes for a specific gene ontology and https://www.biostars.org/p/52101/); following is my code:
ensembl = useMart("ensembl",dataset="hsapiens_gene_ensembl") #retrieve human.ensembl data via biomaRt

go_1 = c("GO:0098830", "GO:0098954", "GO:0099007", "GO:0099067", "GO:0099037", "GO:0098955", "GO:0099592", "GO:0099532") #GO IDs of all terms associated with "presynaptic endosome" (inclusive of child terms)

pre.gene.data <- getBM(attributes=c('hgnc_symbol', 'ensembl_gene_id', 'go_id', 'go_linkage_type'),
                    filters = 'go', values = list(go_1), mart = ensembl)

but my code doesn't filter the genes based on the GO terms. It gives me the following output:
GO:0006886 AP3D1 ENSG00000065000 intracellular protein transport IEA
GO:0016192 AP3D1 ENSG00000065000 vesicle-mediated transport IEA
GO:0030117 AP3D1 ENSG00000065000 membrane coat IEA

... and 49 other entries for the gene "AP3D1"

I'm unable to understand what might be the problem. I replaced the filter "go_id" to "go" because of the updated filters used:listFilters(ensembl).
Please help!

It is filtering by your GO terms, it is only giving you genes associated with those GO terms. However the GO terms in the output are not linked to the GO terms in your input, so it is giving you all the GO terms linked to the genes which are associated with the GO terms you filtered by.
Essentially BioMart is not very clever. It has a table of all the GO terms linked to the genes, it filters the genes to only include those linked to the GO terms you've given it. Then it takes that list of genes and queries that same table completely separately, getting the GO terms for each gene.


Issues installing Dendropy I want to perform phylogenetic analysis using Dendropy. I have installed dendropy using command sudo pip install -U dendropy but when I import it in python3.5.1 it gives followinng error:

 > python3
   Python 3.5.1 |Anaconda 2.5.0 (64-bit)| (default, Dec  7 2015, 11:16:01)
   [GCC 4.4.7 20120313 (Red Hat 4.4.7-1)] on linux
   Type "help", "copyright", "credits" or "license" for more information.
 >>> import dendropy
 Traceback (most recent call last):
 File "<stdin>", line 1, in <module>
 ImportError: No module named 'dendropy'

I again installed dendropy using the above given command as follows:

sudo pip3 install -U dendropy
Requirement already up-to-date: dendropy in /usr/local/lib/python3.5/dist-packages (4.4.0)
Requirement not upgraded as not directly required: setuptools in /usr/local/lib/python3.5/dist-packages (from dendropy) (39.2.0)

After doing this again the error remains the same and I am still not able to import dendropy, any solution to this problem?

You should never install python packages as root or with sudo.

You're already using conda, so conda install -c bioconda -c conda-forge dendropy.


Gap-fill assembly with PacBio reads I would like to gap-fill (or correct) a mammalian-sized assembly for which we have BACs that have been sequenced with PacBio CLR data. I have seen there are pbm22/gcpp tools available for correcting assemblies, but looking at how they've been used, they seem to correct existing assemblies built with FALCON or similar.
Can pbmm2 and gcpp be used with a reference constructed without pacbio data? E.g. given a mammalian-sized reference and a .bam file of PacBio reads for a BAC, what's the best way to gap-fill the current reference? Thx

My understanding is that PBJelly is the tool intended to solve this problem.
This person seems to have used it starting with Illumina scaffolds. I have not personally used it but have heard that it works pretty well.
Edit: note that the PacBio ecosystem is not super easy to use from scratch. However, they have made pretty good efforts to facilitate setup using conda, as for instance documented here. This person shows a conda recipe for pbjelly.


Tools for simulating Oxford Nanopore reads Are there any free open source software tools available for simulating Oxford Nanopore reads?

Simulators designed specifically for Oxford Nanopore:

NanoSim
NanoSim-H
SiLiCO
ReadSim
DeepSimulator

General long read simulators:

Loresim
Loresim 2
FASTQsim
LongISLND

For an exhaustive list of existing read simulators, see page 15 of my thesis, Novel computational techniques for mapping and
classifying Next-Generation Sequencing data.


Help with the definitions of fields in a VCF output by Strelka In Strelka vcf in INFO column we have these

Depth="The number of reads covering the variant position including duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds."

ReadCount="The number of reads covering the variant position excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">
ReadCountControl="The number of reads covering the variant position in the control excluding duplicates, supplementary records and reads that fall below minimum base and mapping quality thresholds.">

And in FORMAT column we have

AU="Number of 'A' alleles used in tiers 1,2">
CU="Number of 'C' alleles used in tiers 1,2">
DP="Read depth for tier1 (used+filtered)">
GU="Number of 'G' alleles used in tiers 1,2">
TU="Number of 'T' alleles used in tiers 1,2">

What is the difference of ReadCount in INFO column and DP in FORMAT column?

What is the difference of number of allele versus number of reads covering the variant position ?

Please help for clarification if you have any intuition on this definitions

Thanks

What is the difference of ReadCount in INFO column and DP in FORMAT column?

I'm not super familiar with Strelka but I believe it follows the VCF spec, so DP in the FORMAT column should give the depth of coverage (ie read count) per sample. ReadCount in INFO is not in the VCF spec but based on the definition it sounds like the depth of coverage aggregated over all samples.

  What is the difference of number of allele versus number of reads covering the variant position ?

I believe the "number of alleles used" refers to counts of each base that are used in the indel error modeling process as described in the preprint here:

https://www.biorxiv.org/content/10.1101/192872v2.full

So it sounds like subsets of the number of reads covering the variant position.


Store contigs and genome substrings based on mismatches I have a string, a list of substrings, and a list of positions where substrings are found in string.

string = 'ATGGTACACGCTACGAGCTAG'
substrings = ['TACA', 'CTAGAAAAG']
positions = ['11, 4', '']

Substring 'TACA' appears at positions 11 and 4 in string (i.e. '11, 4').

Substring 'CTAGAAAAG' does not appear in string (i.e. '').

I need to remove 'CTAGAAAAG' from substrings and I need to remove '' from positions, while creating a new list with 'CTAGAAAAG'. I need your help here to create code which will give the following results:

new_substrings = ['TACA']
new_positions = ['11, 4']
new_list = ['CTAGAAAAG']
string = 'ATGGTACACGCTACGAGCTAG' (obtained from above)

Once this is complete, the given code:

new_positions = [map(int, x.split(',')) for x in new_positions]
seq_pos = dict(zip(new_substrings,new_positions))
seq_ref = {k: [string[start:start+length] for start in starts] for k, starts, length in zip(new_substrings,new_positions,[len(x) for x in new_substrings])}

Gives:

seq_pos = {'TACA': [11, 4]}
seq_ref = {'TACA': ['TACG', 'TACA']}

Then, from seq_ref, I need to create a list of integers which represents the number of mismatches between the key: 'TACA' and both of its values: 'TACG' and 'TACA'.

mismatches = [[1, 0]]

The code to find the # of mismatches between two strings is provided:

def HammingDistance(p, q):
    dist = 0
    for i in range(len(p)):
        if p[i] != q[i]:
            dist += 1
    return dist

Finally, I need to save my data in a CSV file (i.e. Excel or Numbers).

new_positions = ['11, 4'] (obtained above)
new_substrings = ['TACA'] (obtained above)
actual_substrings = [['TACG', 'TACA']] (values of seq_ref)
mismatches = [[1, 0]] (obtained above)

Info needs to be sorted based on # of mismatches.

I understand this is a multi-step problem but any help would be appreciated.

For the first question try something like:

new_substrings = []
new_positions = []
new_list = []
for pos, sub in zip(positions, substrings):
   if not pos:
       new_list.append(sub)
   else:
       new_substrings.append(sub)
       new_positions.append(pos)

this gives:

>>> new_substrings
['TACA']
>>> new_positions
['11, 4']
>>> new_list
['CTAGAAAAG']

For the hamming distance you just need to iterate over the values in seq_ref and append the results of the hamming distance function in mismatches.  Based on the given code and the example above you should be able to create that block.  Below is an untested quick example to get you on track.

mismatches = []
for seq in seq_ref.keys():
   for value in seq_ref[seq]:
       mismatches.append(HammingDistance(seq, value))

As far as the CSV and sorting is concerned I suggest that you look into pandas or numpy.  Both of those packages provided functions to help with this.


How to interpret track height in Chipseq? I recently received some .bigwig files from a chipseq experiment.
I loaded the files into IGV.
The interpretation is rather intuitive.
My question is the following:
What exactly do the track heights represent?
For IgG control, the track height is low ranging from 0~5.
For positive controls, such as H3K4me3, I see clear peaks with height ~100.
Do track heights simply mean the number of reads mapped to that genomic interval during peak calling?
Does it literally mean there are ~100 reads that aligned to that region in the H3K4me3 .bam file?
Thanks so much.

Depends how the files were made. In the simplest case yes, the height represents the pipeup of reads from the BAM file that was used. For a direct (visual) comparison you have to normalize the files though as otherwise sequencing depth confounds the height of the peaks, like if file A has ten times more reads and the peak is 10 times higher than in B then this is probably only due to differences in total read number (=sequencing depth). I wrote down my preferred method to do that in this biostars post: https://www.biostars.org/p/413626/#414440
Does this make sense to you?


Empty .result in MirDeep Star I am using MirDeep star tool for RNA sequence alignment. I have tried both GUI version and command line version of this tool. But every time, I am getting empty . result file and .known_mir file. This tool is neither giving any error nor any warning with my data. I have also tried the demo file given with this and this tool is working fine with this (i.e. giving some results in .result and known_mir file).

Can anyone suggest me what could be the possible reason for this. Why I am getting empty file.

I have found reason for not getting any result. I was using default adapter given by mirdeep* which gives me zero alignment with my data. This error pointed out by author of mirdeep* ( thanks to Jiyuan and colleen). After giving right adapter, I got the right result.


Minimal working example for phylogenetic tree construction in R using an alignment of DNA strings and a vector of names? I am looking at how to build a phylogenetic tree in R from aligned genetic information (dna sequences). There are multiple tutorials on how to go about this task eg tutorial1, tutorial2 ("ape" docs), tutorial3, and there are many more. There appears to be quite a few steps with dependencies from multiple packages: eg. "ape", "phangorn", "phytools", "geiger". I would like to learn a minimal working example (mwe) for constructing a phylogenetic tree from a basic data-structure holding dna-strings.
Eg. I have a vector of nucleotide strings, and a separate vector for the 'names' of each string in the vector, how can I construct a phylogenetic tree from this information? (not a fasta file)

A key difference from many examples provided is that they assume that the data provided is in the form of sequences rather than an already aligned sequence alignment.
library("Biostrings")
library("seqinr")
library("ape")
library(phylogram)
library("dendextend")

#read in the data from a format with the aligned sequences (eg msf)
msf.res <- read.alignment(file = "geneticAlignment.msf", format = "msf")
#dist.alignment and not dna.dist commonly used in tutorials
msf.res.dist.alignment = dist.alignment(msf.res, matrix = "identity")
#perform neighbor joining upon the distance set
msf.res.dist.alignment.nj = nj(msf.res.dist.alignment)
#then plot the joined sequences
plot(msf.res.dist.alignment.nj, main = "Neighbor Joining")


How can I predict if the binding site of a DNA-binding protein is located on the promoter sequence of a gene? Is there any computational tool or any other method to predict if the binding site of DNA binding protein is actually located on a promoter sequence?

You could run FIMO on your entire genome for TFs (transcription factors; DNA-binding proteins) of interest, which gives you predicted binding sites: genomic intervals where those TFs bind, with some statistical significance.

https://bioinformatics.stackexchange.com/a/2491/776

You can then map or intersect promoter regions with those TF binding sites with BEDOPS bedmap or bedops, respectively:

https://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html

https://bedops.readthedocs.io/en/latest/content/reference/set-operations/bedops.html

You can define your promoter regions however you like, e.g., from the TSS to 2.5kb upstream of the TSS. It's up to you and your experiment.

The use of bedmap is recommended to assign TFs to a promoter of interest. This tool will return both the promoter region and any TFs that overlap it by the specified overlap criteria.


Why this ballgown error (first column of pData does not match the names of the folders) comes out I am trying to create a ballgown object, and I would like to set a phenodata information.
For that I do

phenodata <- read.csv('l002.csv')

and my l002.csv file contains the following info

id,state
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube

A ballgown folder contains folders that have exactly the same names, e.g. RNA-CCC_8N_S5_L002, RNA-CCC_8T1_S6_L002, etc...
Then, I use this piece of code for bg object creation

bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata )

but every time I run it, I have the same error every time:

Error in ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata) :
  first column of pData does not match the names of the folders containing the ballgown data.

What can be a problem since first column of pData matches exactly the names of the folders containing the ballgown data..

If I get rid of phenodata ballgown object creation works, but I need this phenodata for further manipulations.

Well, it looks like if you're using Windows for managing files - that's not a good idea - it sorts names of your ballgown folders in alphabetical order while you need it differently.

That is why the one should use Linux - when you're ls -ltr your ballgown folder in Linux it gives you the "right" sequence of folder names to be put inside a ballgown object using pData.

This is what helped me to get rid of this error:

In Linux enter your ballgown folder and do ls -ltr | awk '{print$9}' and you will have the right sequence that you have to put inside you .csv file using Notepad++, indicating what's what. Below's my example compared to the above mentioned initial example:

id,state
RNA-CCC_12N_S7_L002,control
RNA-CCC_12T_S8_L002,tube
RNA-CCC_13N_S105_L002,control
RNA-CCC_13T_S106_L002,tube
RNA-CCC_23N_S107_L002,control
RNA-CCC_23T_S108_L002,tube
RNA-CCC_8N_S5_L002,control
RNA-CCC_8T1_S6_L002,tube

I don't have my error anymore and after

phenodata <- read.csv('my_phenodata_file.csv')
bg = ballgown(dataDir = "ballgown", samplePattern = "RNA-CCC_", pData = phenodata)

my results are nicely wrapped

Sat May 02 22:11:57 2020
Sat May 02 22:11:57 2020: Reading linking tables
Sat May 02 22:11:59 2020: Reading intron data files
Sat May 02 22:12:04 2020: Merging intron data
Sat May 02 22:12:06 2020: Reading exon data files
Sat May 02 22:12:14 2020: Merging exon data
Sat May 02 22:12:16 2020: Reading transcript data files
Sat May 02 22:12:20 2020: Merging transcript data
Wrapping up the results
Sat May 02 22:12:22 2020

The problem is explained here as well

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

Hope this will help anyone


Is it conflicting: significant difference according to Wald test but not in PERMANOVA and ANOSIM to compare species abundances of metagenomic samples? If there is no significant difference between two groups of metagenomic samples according to "Multivariate differential abundance tests" for example PERMANOVA and ANOSIM, does it makes sense to use "Univariate differential abundance testing" to see, for example, what species have significantly different abundances between groups?

The issue is that I found species that have significantly different abundances among samples of two groups using Wald test in DESeq2, but there is no significant difference between the two groups according to PERMANOVA and ANOSIM.

Yes, it does make sense to perform differential abundance testing after one doesn't find compositional differences with PERMANOVA / ANOSIM.

It seems to me you are drawing an (implicit) analogy to ANOVA and post-hoc tests like the Tukey Honest Significant Difference test. Intuitively, the analogy makes sense, however, these situations are quite different. PERMANOVA / ANOSIM and DESeq2 are operating on different aspects of the data, with different transformations and under different assumptions.

PERMANOVA and ANOSIM do not operate on the raw counts directly, rather, they test differences between groups by calculating some distance or dissimilarity between the samples. One can choose between several distance / dissimilarity measures, and the test result may be impacted by the chosen measure.

DESeq2 (and other packages), on the other hand, test for species abundances under a parametric model, operating on the raw counts directly.

Thus (sadly) one is not a natural, optional follow-up of the other, they are complementary and will tell different stories about the data.


Percentage of each cluster in Seurat I am using Seurat to analyze my single cell data. I have 2 conditions, treated and untreated.  I am trying to create a stacked bar graph in order to show the differences in cell types for each condition but need to collect the percentages of each cluster for the specific cell types. How to I put together a sheet that contains the percentages of each cluster and the make a stacked bar graph of comparison of each condition.

pt <- table(Idents(myobj), myobj$orig.ident)
pt <- as.data.frame(pt)
pt$Var1 <- as.character(pt$Var1)

ggplot(pt, aes(x = Var2, y = Freq, fill = Var1)) +
  theme_bw(base_size = 15) +
  geom_col(position = "fill", width = 0.5) +
  xlab("Sample") +
  ylab("Proportion") +
  scale_fill_manual(values = brewer.pal(12, "Paired")) +
  theme(legend.title = element_blank())

This shall work. myobj is your seurat object. This is my old code so if I remember correct, orig.ident is your treated and untreated conditions. Idents should give you cell clusters.
If this does not work, please let me know


Error: package or namespace load failed for 'ChIPseeker' How do I solve the problem of loading the  bioconductor Chipseeker package? I have installed in the following way and I having an error while loading it and I do not know why.

BiocManager::install("ChIPseeker")
BiocManager::install("TxDb.Hsapiens.UCSC.hg38.knownGene")
BiocManager::install("TxDb.Hsapiens.UCSC.hg19.knownGene")

library(ChIPseeker) # produces and error

Error: package or namespace load failed for 'ChIPseeker':
  'TxDb.Hsapiens.UCSC.hg19.knownGene' is not exported from 'namespace: TxDb.Hsapiens.UCSC.hg19.knownGene'

Thank you in advance for any help on how to go about this.

session info:

R version 3.5.1 (2018-07-02)
Platform: x86_64-w64-mingw32/x64 (64-bit)
Running under: Windows 7 x64 (build 7601) Service Pack 1

Matrix products: default

locale:
[1] LC_COLLATE=German_Germany.1252  LC_CTYPE=German_Germany.1252    LC_MONETARY=German_Germany.1252
[4] LC_NUMERIC=C                    LC_TIME=German_Germany.1252

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

other attached packages:
 [1] BiocInstaller_1.32.1                    TxDb.Hsapiens.UCSC.hg19.knownGene_3.2.2
 [3] TxDb.Hsapiens.UCSC.hg38.knownGene_3.4.0 GenomicFeatures_1.34.8
 [5] AnnotationDbi_1.44.0                    GenomicRanges_1.34.0
 [7] GenomeInfoDb_1.18.2                     IRanges_2.16.0
 [9] S4Vectors_0.20.1                        cellrangerRkit_1.1.0
[11] Rmisc_1.5                               plyr_1.8.4
[13] lattice_0.20-35                         bit64_0.9-7
[15] bit_1.1-14                              ggplot2_3.1.1
[17] RColorBrewer_1.1-2                      Biobase_2.42.0
[19] BiocGenerics_0.28.0                     Matrix_1.2-14

loaded via a namespace (and not attached):
Error in x[["Version"]] : subscript out of bounds
In addition: Warning message:
In FUN(X[[i]], ...) : DESCRIPTION von Paket 'digest' fehlt oder zerstört

Simple, use the latest version of R and Bioconductor, now on version 3.6.1 and Bioconductor 3.9. Those are the versions that are supported. If you can't upgrade R, check that BiocManager::valid() is true, you might have ended in some inconsistency between Bioconductor versions.

Looking at your sessionInfo, you seem to have some problems with the installed libraries, check and correct those too if possible


Convert rs ID of one hg build to rs IDs of another build I have a list of dbSNP rsIDs for GRCh37 and I want to convert them to the equivalent IDs in GRCh38. This is using the most recent dbSNP build (150 as of the time of this post). Is there any ID mapping available? If not, are there any tools I can use?

You can assume that the overwhelming majority of rsIDs are the same between GRCh37 and GRCh38 (they're semi-stable IDs). There are, however, a number of rsIDs that are present only in GRCh37, which you can find here. Note that the format of this file is a bit strange, it's chromosome|position|ID|weight, where position is sometimes empty, weight is typically 1, and ID is the rsID without the rs prefix.

Also, some rsIDs have a different strand. This is less likely to be an issue for you, but it's there in case you need it.


pymol script to only select 5 chains within a distance from a reference and save the selection I am generating 30 asymmetric units of a protein structure. Now, I want to select different combinations of 5 units which are 30 Angstroms away from the reference and save it as a different file. Here is the command I could work out till now :

fetch 3q2v
symexp sym,3q2v,(3q2v),30
select chains, bychain all within 30 of 3q2v

How to loop in and select only 5 chains at a time obtained from the setup obtained by the above code, and save it in separate files?

Firstly, I would preface that the PyMOL command line in the GUI is great, but you can run PyMOL as a regular Python library (calling the appropriate python). This makes life so much simpler for more complex manipulations. You can do really complex manipulations with it. I cannot suggest it enough.
The API version of the commands works only marginally differently, but you can do API calls in the GUI cmd.) and PyMOLWiki documentation is only slightly weaker for them (help() fills the gaps).

But if you really want to use the GUI route...

Let's get the names of the models in the selection

onames = cmd.get_object_list('chains')
print onames

Now, you want to save subsets of 5 chains (technically, objects or models). That is, you do not care about the order they are in and you don't want the same one twice (actually selections with 4 or less), so you want all combinations. Make sure to do the following
    print len(onames)

And google that number as <n> choose 5. That is how many files you will get. I am pretty sure most of these sets will be meaningful and will be useless.
In a script a loop is easy, while in the gui you have to keep it simple or use a list coprehension hackishly.

 from itertools import combinations
 for s in combinations(onames, 5): cmd.save('-'.join(s)+'.pdb', ' or '.join(s))

Obviously, you could alternatively use an enumerate method around the iterator and get an index for the file. But do note f-strings will not work.

 for i, s in enumerate(combinations(onames, 5)): cmd.save(str(i)+'.pdb', ' or '.join(s))


Interpreting a kernel density function plot of gene expression I need help interpreting a kernal density plot and a histogram of this data. I have provided some context below let me know if you need more information

A co-worker that I am collaborating with has a list of genes ranked by expression fold change.  He says most of the genes in one dataset are noise with respect to his analysis so he is trying to filter this out.  He has sent me a document of his methods for generating a new list of genes which he think are relevant.  he has run a permutation test on the ranked genes I have provided a quote here:

"To determine which genes contributed more strongly in real data, as compared by permuted data, the rankings achieved for each gene in real data were compared to the shuffled rankings. For each gene, and each permutation type, a p-value was calculated by counting how many times that gene was ranked higher in shuffled data than in real data, and dividing it by the number of comparisons (16 * 256 = 4096). To achieve p<0.05 a gene’s real ranking would have to be higher than 95% of its shuffled rankings."

Following this he provides kernel density plots of the genes under different rankings.

"Below are kernel density estimation (kde) plots, overlaid on histograms, of the rankings of two genes randomly selected from the few genes with p<0.05 vs all three permutations."

"Below are kde plots of the rankings of two genes randomly selected from the many genes with p>0.05 vs all three permutations."

What am I looking at?  What does this tell me about the genes ranking statistics.

There are two topics here to be discussed. Gene Set Enrichment Analysis and a kernel Density plot.
GSEA
The Wikipedia page about Gene Set Enrichment Analysis covers well the principle and the tools that can be used to calculate it. To cite it:

The general steps are summarized below:
Calculate the enrichment score (ES) that represents the amount to which the genes in the set are over-represented at either the top or bottom of the list. This score is a Kolmogorov–Smirnov-like statistic [a ranking based statistic].
Estimate the statistical significance of the ES. This calculation is done by a phenotypic-based permutation test in order to produce a null distribution for the ES. The P value is determined by comparison to the null distribution.
Calculating significance this way tests for the dependence of the gene set on the diagnostic/phenotypic labels

Basically, you want to test if the ranking is random or not. Therefore you make a random dataset and test whether your set is like that. One assumption underlying this is that the output of the gene whatever it may be is not normally distributed (parametric, follows a bell curve), hence the use of ranking, which removes this (you lose signal but you lose biases). I am assuming the non-parametric nature of your ouput values was tested or is a known non-parametric value, such as RPKM or log2 fold change (by virtue of being a ratio). In other words, the magnitude difference between your samples is lost, but their order is not.
There are several tools, some online, that do this enrichment testing for you depending on the task, some even work with gene ontologies.
KDE
A kernel density plot is a like a histogram, but smoothed, albeit not in a moving average way. In a histogram you divide your x-axis into bins: it is discreet and you get an integer count per bin. In a kernel density plot the data is fitted to "probability density function", an equation which given x will give y, where the integral of the curve is one, hence why y is density not counts. The x axis will be your bins-like value, in this case rank. The binning window metric for a KS function is called bandwidth.
See Matlab documentation for a fairly non-opaque explanation.
However, if you want to compare histograms, I strongly suggest violin plots, not because they are strongly in vogue, but because you can compare histograms without worrying about overlaying them.


Gene function annotation - bacterial genome I'm trying to annotate a genome to find all genes with a specific function. I have a FASTA and the read FASTQs - I'd like to assign the functional group of the identified proteins (e.g. Kegg orthology) automatically.

For more context, I have whole genome sequences of bacteria - I'm trying to find all genes involved in membrane transport. I'd either like a tool to retrieve gene ontologies from a list of gene names (input: genA output: Membrane Transport), or a tool that outputs the function as part of the annotation pipeline (e.g. genbank output with GO function defined).

Or posed in a different way, how would you identify all genes with a specific function from prokaryotic whole genome sequence data?

Are there any scripts that can do this?

The PATRICBrc.org is a free resource you can annotate your genome and compare it with publicly available genomes.

Refer this tutorial for more detail.

Once you have annotated your genome in PATRIC, you will have PATRIC Local families (k-mer based protein family), which will help you identify all genes with a specific function.


Why do we scale the rate matrix of a substitution model to make the average rate of substitution equal 1? The matrix of transition probabilities for a substitution model Q over time t is found as follows: $e^{Qt}$. Since Q and t only show up as a product, we cannot distinguish time and rate of change. So we multiple Q by a constant so that the average rate of substitution is 1. Time t is then measured by the expected number of substitutions per site.

But why does forcing the average rate of substitution to equal 1 have this effect?

@Michael's comment is correct but possibly a little difficult to understand, so I'm going to rephrase it here with an example (please let me know if I'm misrepresenting):

This matrix (pretending for a moment that we care only about purine/pyrimidine transitions and thus only two states) is directly interpretable in terms of the number of transitions per unit time:
$$ Q = \begin{matrix}
0 & 1.5 \\
0.5 & 0 \\
\end{matrix}
$$

Whereas this matrix is harder to think about because it is talking about the number of transitions per (I think?) $\frac{1}{0.37}$ units of time:
$$ Q = \begin{matrix}
0 & 0.37 \\
0.37 & 0 \\
\end{matrix}
$$

You will have to rescale any result anyways to talk about it in terms of "per unit time", so you might as well do it in the matrix itself so that at each point we are dealing with a more or less intuitively comprehensible numbers, given that the scaling doesn't actually change anything because time is dimensionless as you point out.

Of course, if you have a time-calibrated phylogeny then time is not dimensionless, so it might make sense to use a different scaling, such that the average rate in $Q$ is instead "per year" or "per million years" or whatever works within the allowances of floating-point precision and brains.

So it's not that rescaling has any particular effect on the math, but rather that rescaling makes it easier for human brains to think about, because we like to represent time in units.


What are the differences between GWAS and QTL mapping? What are the differences between the two methods?

What advantages does one have over the other, and what are their limitations relative to one another?

The basic difference between GWAS and QTL mapping is that GWAS studies the association between alleles and and a binary trait, such as being a sufferer of a disease, while QTL analysis deals with the contribution of a locus to variation in continuous trait like height.

Thus GWAS normally deals with frequencies (i.e. is the frequency disease probands greater in those with the minor allele compared to the major allele at locus X), whereas QTLs deal in correlations (Is the number of minor alleles at locus Y correlated with height).

Methodologically the primary difference would be the distribution that links the allele to the outcome (disese/continuous phenotype). GWAS is going to use Chi-squared statistics or logistic regression, QTL with more likely use straight up normal distribution based linear regression.


What is the best way to drop alignment artifacts in sequence objects in biopython? I have a bunch of sequences created with Roary which I want to preprocess. When I read the sequences (genomics data) I want to translate them to protein sequences using:
import pandas as pd
from Bio.Seq import Seq
from Bio import SeqIO

def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.translate())} for rec in fasta_sequences])
    return sequences

The problem is that .translate() fails with:

TranslationError: Codon 'T--' is invalid

Is there a build-in function to ignore the dashes?

I think the .ungap() is what I was looking for.
def read_fasta_file(fn):
    with open(fn,'r') as file:
        fasta_sequences = SeqIO.parse(file, 'fasta')
        sequences = pd.DataFrame([{
            'filename': basename(fn),
            'id': rec.id,
            'seq': str(rec.seq.ungap().translate())} for rec in fasta_sequences])
    return sequences


Obtaining uniquely mapped reads from BWA mem alignment This is based on a question from betsy.s.collins on BioStars. The original post can be found here.

Does anyone have any suggestions for other tags or filtering steps on BWA-generated BAM files that can be used so reads only map to one location? One example application would be to find seeds for the TULIP assembler/scaffolder, which works best for reads that map to unique genomic locations.

The concept of "uniquely mapped reads" is a loaded term, and most sources suggest filtering by MAPQ should do the trick. However, this approach doesn't seem to work when using BWA as a read mapper.

To exclude all possible multi-mapped reads from a BWA-mapped BAM file, it looks like you need to use grep on the uncompressed SAM fields:

samtools view -h mapped.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > unique_mapped.bam

Explanation follows...

I'm going to assume a situation in which a bioinformatician is presented with a mapped BAM file produced by BWA, and has no way of getting the original reads. One high-effort solution would be to extract the mapped reads from the BAM file and re-map with a different mapper that uses the MAPQ score to indicate multiple mappings.

... but what if that were not possible?

My understanding of BWA's output is that if a read maps perfectly to multiple genomic locations, it will be given a high mapping quality (MAPQ) score for both locations. Many people expect that a read that maps to at least two locations can have (at best) a 50% probability of mapping to one of those locations (i.e. MAPQ = 3). Because BWA doesn't do this, it makes it difficult to filter out multiply-mapped reads from BWA results using the MAPQ filter that works for other aligners; this is likely to be why the current answer on Biostars [samtools view -bq 1] probably won't work.

Here is an example line from a BWA mem alignment that I've just made. These are Illumina reads mapped to a parasite genome that has a lot of repetitive sequence:

ERR063640.7     16      tig00019544     79974   21      21M2I56M1I20M   *       0       0       TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13      AS:i:77 XS:i:67 XA:Z:tig00019544,-78808,21M2I56M1I20M,6;tig00019544,-84624,79M1I20M,6;tig00019544,-79312,33M4I42M1I20M,8;

BWA mem has found that this particular read, ERR063460.7, maps to at least three different locations: tig00019544, tig00019544, and tig00019544. Note that the MAPQ for this read is 21, so even though the read maps to multiple locations, MAPQ can't be used to determine that.

However, the alternative locations are shown by the presence of the XA tag in the custom fields section of the SAM output. Perhaps just filtering on lines that contain the XA tag will be able to exclude multiply-mapped reads. The samtools view man page suggests that -x will filter out a particular tag:

$ samtools view -x XA output.bam | grep '^ERR063640\.7[[:space:]]'
ERR063640.7 16  tig00019544 79974   21  21M2I56M1I20M   *   0   0   TATCACATATCATCCGACTCAGCTCGACGAGTACAATGCTAATTTAACACTTAGAATGCCCGGCAATGAAATTCGTTTTCCGTCAATTCTTGAAAATTTC    <AABBEGABFJKKKIM7GHKKJK>JLKLDGMHLIMIHHCGIJKKLJKLNJGLLLKLILKLMFNDLKGHJEKMKKMIJHGLOJLLLKIJLKKJEJLIGG>D    NM:i:4  MD:Z:83A13  AS:i:77 XS:i:67

... so it filtered out the tag (i.e. the tag no longer exists in the SAM output), but not the read. There are no useful bits in the FLAG field to indicate multiple genomic mappings (which I know can be filtered to exclude the read as well), so I have to resort to other measures.

In this particular case, I can use grep -v on the uncompressed SAM output to exclude alignment lines that have the XA tag (and re-compress to BAM afterwards, just to be tidy):

$ samtools view -h output.bam | grep -v 'XA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | grep '^ERR063640\.7[[:space:]]'
<no output>

Hurray! reads filtered. As a little aside, this grep search has a fairly substantial computational load: it's looking for some string with the text XA:Z: somewhere in the line, and doesn't actually capture every situation. Some masochistic person might come along at a later date and decide that they're going to call all their reads HAXXA:Z:AWESOME!:<readNumber>, in which case a tweak to this grep search would be needed to make sure that there's a space (or more specifically a tab character) prior to the XA:Z:.

Now I do a check for any duplicated read names, just to be sure:

$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
ERR063640.1194
ERR063640.1429
ERR063640.1761
ERR063640.2336
ERR063640.2825
ERR063640.3458
ERR063640.4421
ERR063640.4474
ERR063640.4888
ERR063640.49
ERR063640.4974
ERR063640.5070
ERR063640.5130
ERR063640.5300
ERR063640.5868
ERR063640.6116
ERR063640.6198
ERR063640.6468
ERR063640.6717
ERR063640.6797
ERR063640.7322
ERR063640.750
ERR063640.7570
ERR063640.7900
ERR063640.8115
ERR063640.8405
ERR063640.911
ERR063640.9206
ERR063640.9765
ERR063640.9986

Oh... damn. I wonder what they are:

$ samtools view output_filtered.bam | grep '^ERR063640.3458[[:space:]]'
ERR063640.3458  16  tig00002961 5402    60  58S38M  *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATATCCGTATCGCGAAGATCAGGACTTACTCCGCAGAAGAA    DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJDHLNKCK;KMKGMFKJILIIIMKI9JLKKHEJFII?CC    NM:i:0  MD:Z:38 AS:i:38 XS:i:0  SA:Z:tig00002377,202353,-,14M3I5M1I35M38S,19,5;
ERR063640.3458  2064    tig00002377 202353  19  14M3I5M1I35M38H *   0   0   AGGTACCATTCGATAGAGGGAGAAAGGCACTACTAAAGATTTTGCCACATTTGCTATA  DD6HFFJBKFH=KDILKLGLJEKLKGFJIH8IKHLLMJEK:L:HBGJIHJKFLLKIHJ  NM:i:5  MD:Z:5G48   AS:i:35 XS:i:27 SA:Z:tig00002961,5402,-,58S38M,60,0;

Aha! Supplemental alignments, which use the official SA tag [other canonical alignments in a chimeric alignment]. These appear to be situations where a single read has been split up, and maps to multiple locations. Note that in this case, both of the alignments still have MAPQ scores of over 3. It sounds like the questioner would also want to get rid of these situations as well. This time, there are standard flag fields as well to deal with these situations (0x800: secondary alignment). Except it's not enough to just filter the supplemental alignment, because both read mappings should be removed, rather than just the one (or ones) that happened to be tagged as secondary.

Luckily, BWA appears to put the SA tag into all reads containing supplementary alignments (if this is not the case, I'm sure someone will correct me on that). So, I add in the SA search as an additional grep filter:

$ samtools view -h output.bam | grep -v -e 'XA:Z:' -e 'SA:Z:' | samtools view -b > output_filtered.bam
$ samtools view output_filtered.bam | awk '{print $1}' | sort | uniq -d
<no output>

Done. Easy peasy! </sarcasm>

... that original "high-effort" solution of using a different aligner doesn't look so bad now.


Error: sort order of reads in BAMs must be the same I am a novice to bioinformatics trying to run a cuffdiff job in cufflinks and keep getting this error message:

Error: sort order of reads in BAMs must be the same

GENERATING SAM FILES:

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_B73_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_Zea_mays.AGPv4.dna.toplevel -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11-5_FCHHL7JBBXX_L7_DLPRMzP300H3AAAPEI-27_2.fq

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11R-5
module load hisat2
hisat2 -k 1 -p 8 -S 11-5_Mo17_accepted_hits.sam -x /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2//Mol17_Ref_Lai -1 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_1.fq  -2 /fs/project/PMIU0157/Claire_Stoll/RNASeq/11R-5_DLPRMzP300B4AAAPEI-2_2.fq

CONVERTING SAM TO BAM:

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load samtools
samtools view -bS a_B73_11-5_accepted_hits.sam > a_B73_11-5_accepted_hits.bam
samtools sort a_B73_11-5_accepted_hits.bam -o a_B73_11-5_accepted_hits.sorted.bam

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5
module load samtools
samtools view -bS a_Mo17_11-5_accepted_hits.sam > a_Mo17_11-5_accepted_hits.bam
samtools sort a_Mo17_11-5_accepted_hits.bam -o a_Mo17_11-5_accepted_hits.sorted.bam

CUFFDIFF JOB

cd /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5
module load cufflinks/2.2.1
cuffdiff -b /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5 -u -N -p 8 /fs/project/PMIU0157/Claire_Stoll/Zea_mays.AGPv4.36.chr.gff3 /fs/project/PMIU0157/Claire_Stoll/RNASeq/B73_hisat2/B73_hisat2_11-5/a_B73_11-5_accepted_hits.sorted.bam /fs/project/PMIU0157/Claire_Stoll/RNASeq/Mo17_hisat2/Mo17_hisat2_11-5/a_Mo17_11-5_accepted_hits.sorted.bam

I think cuffdiff is not quite giving you the right error...it looks like the crux of your problem is you aligned each sample to a different reference.


the model matrix is not full rank : This is a classic question which a biologist face without clear understanding of the model design Saw this answer biostar. Tried to make my metadata as such but still the
"Error in checkFullRank(modelMatrix) : "

This is my coldata
dput(coldata)
structure(list(Group = structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L,
2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L,
3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 1L,
1L, 1L, 1L, 1L, 1L, 1L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 6L, 6L,
6L, 6L, 6L, 6L), .Label = c("HSC", "Blast", "CMP", "GMP", "LSC",
"Mono"), class = "factor"), Mutation = structure(c(10L, 4L, 3L,
5L, 2L, 3L, 9L, 11L, 1L, 1L, 7L, 7L, 10L, 7L, 8L, 2L, 10L, 10L,
10L, 10L, 2L, 2L, 7L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L,
6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 6L, 2L, 2L, 5L, 2L,
9L, 1L, 1L, 8L, 6L, 6L, 6L, 6L, 6L, 6L), .Label = c("DNMT3A-R882H",
"FLT3-ITD", "IDH2-R140Q", "KRAS-G13D", "MAU2-Q98H", "No", "NPM1-L287ins(TCTG)",
"NRAS-G13D", "STAG2-R614", "TET2-E1357", "TET2-R550"), class = "factor"),
    Satus = structure(c(2L, 5L, 4L, 6L, 3L, 4L, 9L, 10L, 1L,
    1L, 7L, 7L, 2L, 7L, 8L, 3L, 2L, 2L, 2L, 2L, 3L, 3L, 7L, 16L,
    16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L, 16L,
    16L, 16L, 17L, 17L, 17L, 17L, 17L, 17L, 17L, 12L, 12L, 13L,
    12L, 9L, 11L, 11L, 14L, 15L, 15L, 15L, 15L, 15L, 15L), .Label = c("Blast-DNMT3A-R882H",
    "Blast-E1357", "Blast-FLT-ITD", "Blast-IDH2-R140Q", "Blast-KRAS-G13D",
    "Blast-MAU2-Q98H", "Blast-NPM1-L287ins(TCTG)", "Blast-NRAS-G13D",
    "Blast-STAG2-R614", "Blast-TET2-R550", "LSC-FLT-DNMT3A-R882H",
    "LSC-FLT-ITD", "LSC-MAU2-Q98H", "LSC-NRAS-G13D", "Mature-Normal",
    "Progenitor-Normal", "Stem-Normal"), class = "factor")), row.names = c("Blast1",
"Blast10", "Blast11", "Blast12", "Blast13", "Blast14", "Blast15",
"Blast16", "Blast17", "Blast18", "Blast19", "Blast20", "Blast2",
"Blast21", "Blast22", "Blast23", "Blast3", "Blast4", "Blast5",
"Blast6", "Blast7", "Blast8", "Blast9", "CMP1", "CMP2", "CMP3",
"CMP4", "CMP5", "CMP6", "CMP7", "CMP8", "GMP1", "GMP2", "GMP3",
"GMP4", "GMP5", "GMP6", "GMP7", "HSC1", "HSC2", "HSC3", "HSC4",
"HSC5", "HSC6", "HSC7", "LSC1", "LSC2", "LSC3", "LSC4", "LSC5",
"LSC6", "LSC7", "LSC8", "Mono1", "Mono2", "Mono3", "Mono4", "Mono5",
"Mono6"), class = "data.frame")

dds <- DESeqDataSetFromMatrix(countData=df2, colData=coldata, design= ~ Group + Mutation)

This is the design I'm trying to incorporate. What i get is my design is nested in one of the column by looking at numerous example.
What i would like to do is I'm trying to run a LRT test on this dataset which works fine if i give "Group" as my design and then i run this
dds_lrt <- DESeq(dds, test="LRT", reduced = ~ 1)

which gives me diffrenitaly expressed across multiple groups. Now there are two question

Is it possible to include interaction such as in my case Group+Mutation and do the lrt test?
Or the above is conceptually wrong to do ?

> table(doldat$Group, coldat$Mutation)

        DNMT3A-R882H FLT3-ITD IDH2-R140Q KRAS-G13D MAU2-Q98H No NPM1-L287ins(TCTG) NRAS-G13D STAG2-R614 TET2-E1357 TET2-R550
  HSC              0        0          0         0         0  7                  0         0          0          0         0
  Blast            2        4          2         1         1  0                  4         1          1          6         1
  CMP              0        0          0         0         0  8                  0         0          0          0         0
  GMP              0        0          0         0         0  7                  0         0          0          0         0
  LSC              2        3          0         0         1  0                  0         1          1          0         0
  Mono             0        0          0         0         0  6                  0         0          0          0         0

Look at this table, you see by eye how some mutations are nested with Group. If you want to include it then you will need to subset the experiment into multiple separate runs, depending on the question you want to answer.


Picking number of samples for clustering I have expression data for 100 samples and clinical data for 50 among them. I want to carryout k-means clustering using expression data and then do a survival analysis (Kaplan-Meier analysis) using the labels obtained from k-means.
Can I used the expression data of all 100 samples to carryout clustering or should I use only 50 samples for clustering?

In case there is no systematic difference between the samples with or without the clinical data (i.e. one set is microarray and the other is RNA-seq), I would use all of the samples. And if the cluster labels obtained with "clustering 50 samples" and those obtained with "clustering 100 samples" does not differ too much, then doing either of two would not have much effect anyway. If you would like to compare cluster labels in that way and assuming you would want to test a number of ks, it would be better to compare each combination of ks and number of samples.


Binding columns with matched elements I have a matrix in which sample IDs are in rows and genes are in column like below

> head(data)
                                 ETV7    PLAUR  SLC11A1
TCGA.JY.A6FB.01A.11R.A336.31 1.791167 6.241607 3.709848
TCGA.JY.A938.01A.11R.A37I.31 5.389175 7.048782 3.124435
                                BIRC3 TNFRSF1B
TCGA.JY.A6FB.01A.11R.A336.31 6.726991 5.758026
TCGA.JY.A938.01A.11R.A37I.31 8.097783 6.621900
>

In another matrix I have sample IDs, months to death and status while the order of sample IDs is not the same with my previous matrix like below

Patient_ids OS_time  OS_event
TCGA-IG-A4QS    3.933333333 1
TCGA-JY-A6FB    61.23333333 0
TCGA-JY-A938    35.33333333 0
TCGA-JY-A939    22  0
TCGA-JY-A93C    23.5    0
TCGA-JY-A93E    25.56666667 0
TCGA-L5-A4OJ    21.3    0
TCGA-L5-A4OS    59.4    0
TCGA-L5-A8NM    7.866666667 1
TCGA-R6-A6KZ    5.133333333 1
TCGA-R6-A6XG    38.93333333 0

How I can bind OS_time and OS_event columns from second matrix to the first matrix keeping the right sample order?

  EDITED

> # Reproducible example:
> patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
> ETV7 = c(1.79,5.38,5,6)
> PLAUR = c(6.24,7.05,5,6)
> SLC11A1 = c(3.7,3.12,5,6)
> df = data.frame(ETV7,PLAUR, SLC11A1)
> rownames(df) = patient_id
> # OS data
> Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
> OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
> os_data = data.frame(Patient_ids,OS_time)
> library(rebus.base)
> df$Patient_id = gsub(DOT,"-",rownames(df))
> position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
> position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
> os_data$Position_matched = unlist(position)
> m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
> colnames(m) = colnames(df)
> for(i in 1:length(position))
+ {
+   if(is.na(position[i])){}
+   else{
+     m[i,] = df[position[i],]
+   }
+ }
Error in xj[i] : invalid subscript type 'list'
>

Here a little way to do it. I don't think it is a perfect solution but it should work.
The trick here is to compare patients_id that are wrote with a . or with a -.
So, here is how I do (I modified your example in order to reflect several configurations):

# Reproducible example:
patient_id = c("TCGA.JY.A6FB.01A.11R.A336.31","TCGA.JY.A938.01A.11R.A37I.31","TCGA.L5.A4OS.01A.11R.A212.31","TCGA.AA.A400.01A.11R.A212.31")
ETV7 = c(1.79,5.38,5,6)
PLAUR = c(6.24,7.05,5,6)
SLC11A1 = c(3.7,3.12,5,6)

df = data.frame(ETV7,PLAUR, SLC11A1)
rownames(df) = patient_id

# OS data
Patient_ids = c("TCGA-IG-A4QS","TCGA-JY-A6FB","TCGA-JY-A938","TCGA-JY-A939","TCGA-JY-A93C","TCGA-JY-A93E","TCGA-L5-A4OJ","TCGA-L5-A4OS")
OS_time = c(3.93,61.23,35.33,22,23.5,25.56,21.3,59.4)
os_data = data.frame(Patient_ids,OS_time)

Now, I create a new column in df that will contains patient_id wrote with a -:

library(rebus.base)
df$Patient_id = gsub(DOT,"-",rownames(df))

Now, I'm looking for common pattern between df$Patient_id and os_data$Patient_ids, find their matching row position and extract these lines:

position = sapply(os_data$Patient_ids,function(x){grep(x,df$Patient_id)})
position = lapply(position,function(x){c(x,rep(NA,1-length(x)))})
position = unlist(position)
os_data$Position_matched = position

Now, we are creating a matrix m that we will fill with matching values

m = data.frame(matrix(ncol = length(df),nrow = nrow(os_data)))
colnames(m) = colnames(df)

for(i in 1:length(position))
{
  if(is.na(position[i])){}
  else{
    m[i,] = df[position[i],]
  }
}

And now, we can bind both os_data and m:

New_df = cbind(os_data,m)

And here the final output:

> New_df
   Patient_ids OS_time Position_matched ETV7 PLAUR SLC11A1                   Patient_id
1 TCGA-IG-A4QS    3.93               NA   NA    NA      NA                         <NA>
2 TCGA-JY-A6FB   61.23                1 1.79  6.24    3.70 TCGA-JY-A6FB-01A-11R-A336-31
3 TCGA-JY-A938   35.33                2 5.38  7.05    3.12 TCGA-JY-A938-01A-11R-A37I-31
4 TCGA-JY-A939   22.00               NA   NA    NA      NA                         <NA>
5 TCGA-JY-A93C   23.50               NA   NA    NA      NA                         <NA>
6 TCGA-JY-A93E   25.56               NA   NA    NA      NA                         <NA>
7 TCGA-L5-A4OJ   21.30               NA   NA    NA      NA                         <NA>
8 TCGA-L5-A4OS   59.40                3 5.00  5.00    5.00 TCGA-L5-A4OS-01A-11R-A212-31

Warning: This script should work for some simple case like yours but take extra cautious to the grep output as the short patient_id term of the os_data could eventually generate several hits instead of one.


ROC analysis, CAFA 2 experiment My question is about ROC curves used in CAFA2 experiment. In this paper they used the ROC analysis for the term-centric evaluation. In order to
perform ROC curve analysis we should have a continuous variable and a classifier (categorical) variable. In CAFA2 experiment,
which variables did they consider as a categorical and continuous? Thank you.

The paper doesn't seem to have been written for general readers. Unfortunately I don't understand the equations enough to be able to explain what the categorical and continuous variables are. Just in case it helps, here are the equations from the paper showing the precision/recall calculations:

$$\begin{array}{@{}rcl@{}} \text{pr}(\tau) &=& \frac{1}{m(\tau)}\sum\limits_{i=1}^{m(\tau)} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \wedge f \in T_{i}\right)}{\sum_{f} \unicode{x1D7D9}\left(f \in P_{i}(\tau) \right)},\\ \text{rc}(\tau) &=& \frac{1}{n_{e}}\sum\limits_{i=1}^{n_{e}} \frac{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in {P}_{i}(\tau) \wedge f \in T_{i}\right)}{{\sum\nolimits}_{f} \unicode{x1D7D9}\left(f \in T_{i} \right)}, \\ F_{\max} &=& \max_{\tau} \left\{ \frac{2\cdot \text{pr}(\tau)\cdot \text{rc}(\tau)}{\text{pr}(\tau) + \text{rc}(\tau)} \right\}, \end{array} $$

  where P i (τ) denotes the set of terms that have predicted scores greater than or equal to τ for a protein sequence i, T i denotes the corresponding ground-truth set of terms for that sequence, m(τ) is the number of sequences with at least one predicted score greater than or equal to τ, $\unicode{x1D7D9}\left (\cdot \right)$ is an indicator function, and n e is the number of targets used in a particular mode of evaluation.


Simulating phenotype with the 1000 Genomes Project I'm looking for a way to simulate phenotypes against a real SNP data source, such as the 1000 Genomes. It must be free for commercial purpose (Eg.: MIT license). Any recommendation?
I'm trying to use the GCTA64, but I couldn't get it working. The documentation doesn't help much.
At the end of the day, I want:
1 - Simulate Case/Control and/or quantitative phenotype
2 - Link it with a real SNP dataset (eg: 1000 Genomes)
3 - Conduct GWAS analysis using Plink and/or Hail.

Looking Hail's and Plink's tutorials for GWAS, I realised both use simulated phenotypes from real SNP data sources (1000 Genomes and HapMap, respectively), but how they created the datasets are beyond of the scope of the tutorials, thus not reported.

As mentioned before, I've tried gcta64, but no success. Here's what I've tried:

1 - Downloaded 1000 Genome sample from Plink page:
Entire dataset as a single .tar.gz (1.12 GB) (A2 allele major, not ref, on chr3 before 15 Oct 2017)

2 - Tried to generate the simulate data by:

./gcta64  --bfile 1kGenomesP1/1kg_phase1_all/1kg_phase1_all --simu-qt  --simu-causal-loci causal.snplist  --simu-hsq 0.5 --simu-rep 3  --keep test
.indi.list --out 1kg_phase1_all

Error: Error: --keep test.indi.list not found.

Question: What should be the files:causal.snplist, test.indi.list, teest? I know the property name suggests what should be, but there's no example of file structure nowhere... Any practical example or tutorial?

Btw - I apologise in advance if this is a too trivial question. I'm quite new at it. Appreciate your patience and help :)

I managed to get it done without the gcta64.
I have written a Jupyter Notebook with the steps: https://github.com/bambrozio/bioinformatics/blob/master/utils/1k-genomes-phenotype-simulator.ipynb


How do I find the most similar sequences from a large set of short sequences? I've been asked to evaluate whether the MinION will be sufficient to distinguish between 20bp CRISPR guide RNAs from the GeCKO v2 set. I know that this set has some exactly identical sequences, which wouldn't be distinguishable regardless of the sequencing method used, so I would like to know whether any two non-indentical sequences differ by less than 15% of the 20bp (i.e. only 3 bases different).

There are around 130k 20bp sequences in the dataset, and I'm pretty sure that an all-vs-all approach is not feasible, but perhaps there's some trick that I'm missing.

What I thought of doing was looking at dimer and trimer/codon counts to hunt down sequences that were most similar. I found about 100 sequences that shared the same trimer count signature (but were not identical), and within this group found [only] one pair with 85% identity; all the remainder were 70% or less.

... but I'm having trouble ignoring the doubt I have that these two sequences are the most similar. It might be possible that there are sequences in the dataset that don't have identical trimer count signatures, but differ by less than three bases. I'll give a very obvious example (which isn't in the dataset):

1: AAAAA AAAAA AAAAA AAAAA
2: AAAAC AAAAA AAAAA AAAAA

The trimer counts for these two sequences are different (Sequence 1 -- AAA: 18; Sequence 2 -- AAA: 15; AAC: 1; ACA: 1; CAA: 1), and yet the sequences are only 1 base different (which might be a problem for the MinION). How can I work out if any such sequence pairs exist in the dataset?

130k * 20bp is a small data set. At this scale, SSE2 Smith-Waterman may work well:

git clone https://github.com/attractivechaos/klib
cd klib && gcc -O2 -D_KSW_MAIN ksw.c -o ksw -lz
./ksw -a1 -b1 -q1 -r1 -t14 20bp.fa 20bp.fa > out.tsv

I simulated 130k 20bp reads from E. coli K-12. The last command line takes about ~2 CPU hours based on the partial output.

It is possible to make this much faster by using 5bp seeds. You can also try ssearch and fasta from the fasta aligner package. They are slower, probably because they are not optimized for such input. Computing Hamming distance with SSE2 will be much faster, but is probably of little use as it does not allow gaps. There are also Gene Myers' edit-distance based algorithms, which can be faster than SSE2-SW.

Bowtie2 (as well as most short-read mappers) won't work well at the default setting. Bowtie2 uses long exact seeds. It will miss many 1-mismatch hits, let alone 3-mismatch ones. The number you get from bowtie2 is an underestimate. You might be able to tune bowtie2, but to find 3-mismatch hits, you have to make it a lot slower.


Is it possible to correct bond order and conectivity problems in an sdf file according to mol2 files? I was doing a Virtual Screening Search using Pharmer software. The .sdf output file contains a column in mol format and another column with Root Mean Square Deviations. I can see some bond order and connectivity problems in the mol column. I do have a mol2 file that contains the original molecules and I was wondering if I can use these mol2 files, to fix this. Can rdkit work in this case?

RDKit has a one line command to fix the bond order (docs)
 mol = Chem.MolFromMolFile('wonky.mol')
 template = Chem.MolFromMol2File('template.mol2')
 newMol = AllChem.AssignBondOrdersFromTemplate(template, mol)

The problem is that the two need to be identical. I see a carboxylate in your screenshots, so you most likely did the correct thing and docked with the appropriate protonation for pH 7, i.e. as the conjugate base with a negative formal charge. If one does not then RDKit will fail.
Also, RDKit is "not great" with mol2 files. If you have the original SMILES that is better. If the problem is atom names, using a PDB is a better call. If you want to assign atom names, this blog post of mine may help.


Comparing two multi-fasta files of the same set of proteins with parser - to find and count mutations after treatment My task is to count the mutations  occurred in several proteins after a treatment. The sequences are all present in the two files in the same order. I opened both files with the FASTA parser (SeqIO.parse) in Biopython and I got all the proteins listed (separated before and after the treatment).

How can I zip the parsers together to count the mutations?

How can I count the mutations that occured after the treatment?

from Bio import SeqIO
for normal_samples in SeqIO.parse("/data/statistic/normal_samples", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

for treated_samples in SeqIO.parse("/data/statistic/with_treatment", "fasta"):
    print(normal_samples.id)
    print(repr(normal_samples.seq))
    print(len(normal_samples))

dict_n_t = dict(zip(normal_samples & treated_samples))

from Bio import SeqIO
for normal, treated in zip(SeqIO.parse("/data/statistic/normal_samples", "fasta"),
                           SeqIO.parse("/data/statistic/with_treatment", "fasta")):
    ... do stuff...

That's generally how you zip iterators together in python.


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

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

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

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

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

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


Tool to remove a PCR contamination in NGS data BACKGROUND
I work on NGS data (illumina paired ends reads) coming from a full extract of RNA (metagenomic). We are interested in the viral fraction of this extract.
I observed a contamination with a PCR amplicon. This contamination has been identified as coming from a PCR result that we used to make near the library preparation. So I know exactly the primers used in the PCR.
Here you can see a tablet visualization illustrating the contamination (2 contaminations). I obtained that image with a mapping on a known reference of my virus.
PROBLEM
The problem I have is that these PCR products are exactly the type of sequence I want from my extract (identifying part of the virus). So I have some "false positive" and some misidentification of viruses present in the extract.
I have thought about removing all reads containing primer sequences but it leads to a gap in my data and I can't identify the viruses in my extract anymore.
I am thinking that all of the read pairs that begin before or end after my primers are "true-positives" and I would really like to keep them (even having the primer sequence in them).

Do you know any tool that can help me in that situation? i.e. removing reads containing a sequence but only if it begins with it.
If the tool can also handle the degeneration of my primer it would be perfect.
An example of primer I used is
CNTGGGAGGGCGATCGCAA

Complementary information (asked in comments)
I have around 30 million read pairs (60 millions reads in total). Contamination reads are between 200.000 to a 1 millions reads (depending on sample). "true" reads (not from the contamination) represent around 0.5% of the contamination.

This bash script will get you all of the reads that do not start with your primer sequence. It also handles degeneracy. It requires the tool bioawk.

It works by looking at the first bases in the read and only prints out the headers that do not have a match with the degenerate primer. Then, comm makes sure that the read was retained in both the R1 and R2 filtering. This allows reads to be filtered whether or not the primers start at the beginning of the R1 or R2 read.

seqtk subseq is where you recover the reads that don't start with the primer.

Be warned though, any off-by-one or sequencing errors will allow some of the sequences with primers at the beginning through. Also, this script requires that the headers for the R1 and R2 reads are the same and aren't appended with anything like "\1", "\2", et cetera.

#!/bin/bash
get_primer_free_starts () {
    bioawk -cfastx 'BEGIN{seqs[0] = "CATGGGAGGGCGATCGCAA"; seqs[1] = "CCTGGGAGGGCGATCGCAA"; \
        seqs[2] = "CGTGGGAGGGCGATCGCAA"; seqs[3] = "CTTGGGAGGGCGATCGCAA"} \
        {thisseq = substr($seq, 1, length(seqs[0])); printme=1\
    for (i=0; i<=3; i++) { \
        if (thisseq != seqs[i]){ \
            printme=0 \
        }
    }; if (printme==1){ \
        print($seq)
    }}' "${INPUT} | sort > "${OUTPUT}
}

INPUT="R1.fastq.gz"
OUTPUT="R1_OK.txt"
get_primer_free_starts

INPUT="R2.fastq.gz"
OUTPUT="R2_OK.txt"
get_primer_free_starts

comm -12 R1_OK.txt R2_OK.txt | sort | uniq > passing_reads.txt

seqtk subseq R1.fastq.gz passing_reads.txt > R1.filtered.fq
seqtk subseq R2.fastq.gz passing_reads.txt > R2.filtered.fq


Where can I find a database that has phenotype information together with associated SNPs? I have been searching for a long time and the furthest I got is some database with the functional description for genes. Then, I have to parse these descriptions manually to figure out the association. Moreover, for some species, there are not even any functional descriptions.

I need this information for mice, rice, and Arabidopsis Thaliana. Can anyone help me?

International Mouse Phenotyping Consortium is building a database of phenotypes and knock-outs of mouse. I believe that this database will be fairly complete (20000 knock-outs), but these are knock-outs, not SNPs... There are several mouse GWAS studies, but I am not aware of a database that would pull all the results together.

Arabidopsis big GWAS project is called 1001 genomes project. A lot more information about Arabidopsis from individual experiments is collected and maintained by The Arabidopsis Information Resource, however this resource is subscription based, therefore what you need might not be available for free.

As mentioned, SNPs of rice were identified and linked to different rice groups in The 3,000 rice genomes project. There is also GWAS done on 512 individuals, not sure if it is sufficient for reliable conclusions (even they claim that they found supersigificant links).

Is this what you were searching for?


Counting number of possible alignments between two DNA sequences using python I was looking BioPython to calculate the number of possible alignments between two sequences. Let's say there are two input sequences with the length of m and n. The program will count the possible number of alignments based on the length of two sequences.

Is it possible to do that with BioPython? If not, how can I do that?

A closed-form solution is offered in An exact formula for the number of alignments
between two DNA sequences by Torres, Cabada, and Nieto:

$$f(m,n)=\sum_{k=0}^{min(m,n)}2^{k}\binom{m}{k}\binom{n}{k}$$

If this solution seems reasonable, you could calculate this without BioPython, but with a simple xrange loop, the power operator, and scipy.special.binom:

#!/usr/bin/env python

import sys
import scipy.special

def alignments(m, n):
    if m < 0 or n < 0:
        raise ValueError('m and n should be non-negative')
    s = 0
    b = min(m, n) + 1
    for k in xrange(0, b):
        s += (2 ** k) * scipy.special.binom(m, k) * scipy.special.binom(n, k)
    return s

def main():
    m = 4
    n = 2
    a = alignments(m, n)
    assert (m == 4 and n == 2 and a == 41), "something went wrong!"
    sys.stdout.write("alignments(%d, %d) = %d\n" % (m, n, a))

if __name__ == "__main__":
    main()

I don't know what your m and n are, but keep in mind that this solution grows quickly. Also, I did some digging and scipy approximates the calculation of the binomial coefficient, which is probably going to be an issue for large m and n.

If you're working with large values, you may want to look into using the Python long type and consider the accuracy of the answer. The Stack Overflow link offers alternatives to scipy that generate an exact calculation of the binomial coefficient, at the cost of a longer runtime.


Can I export the SNPs between whole genome alignments using the command line? I am using Mauve to align two whole genomes. I use the following command to get the alignment in xmfa format:

/Applications/Mauve.app/Contents/MacOS/progressiveMauve ref1.fa ref2.fa --output alignment.xmfa

I then open the xmfa file in the Mauve GUI and export the SNPs. Is there a way of doing this using a command? I have tried to follow the procedure described here but the link to snapshot releases of Mauve does not work anymore.

I could not find anything in the Mauve command line help:

/Applications/Mauve.app/Contents/MacOS/progressiveMauve --help

nor on the Mauve website.

I tried the procedure described by Konrad Rudolph. I wrote:

MAUVE_DIR=/Applications/Mauve.app/Contents/Java

in the terminal. Is this correct? Don't I need an export and quotation marks to create a bash variable? I tried adding an export and using quotation marks for the file but the error below remained the same.

followed by:

export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"

followed by:

java org.gel.mauve.analysis.SnpExporter -f alignment.xmfa -o aln.snps

and it seems to work.

It works for me with the download from download page. The download archive contains Mauve.jar so the command line in the answer you link to should work.

However, you might need to include all the jar files (not just Mauve.jar) in the classpath for this to command work, though. For instance, the following works for me:

MAUVE_DIR=‹the-dir-containing-Mauve.jar›
export CLASSPATH="$(find "$MAUVE_DIR" -name \*.jar -print0 | tr '\0' :)$CLASSPATH"
java org.gel.mauve.analysis.SnpExporter -f aln.xmfa -o aln.snps

… there’s probably a smarter way to make this work but I don’t use Java enough to know how.


DESeq2 compare within a condition This might be a really stupid question, but I can't figure out how to do this. I've read through the DESeq2 vignette and manual pages but couldn't find an answer.

I have a bunch of samples split up into different conditions (eg. celltypes and disease state). I would like to make a comparison between the two possible values within one condition, but only for the samples with a specific value of the other condition.

For example, with the following samplesheet:

patient | phenotype | type
--------+-----------+-----
1       | healthy   | A
1       | healthy   | B
1       | sick      | A
1       | sick      | B
2       | healthy   | A
2       | sick      | A
2       | sick      | B

I would like to compare "healthy" vs "sick", but only for type "A".

Currently, I have the following code, but this will also include the type "B" sample:

dds <- DESeqDataSetFromMatrix(countData=counts, colData=design, design = ~ patient + phenotype + type)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("phenotype", "healthy", "sick"))

Any idea how to accomplish this?

Welcome to the world of all possible pairwise comparisons in DESeq2 !

So the easiest way is to create another group in your data frame, first I simulate some sensible counts:

counts = counts(makeExampleDESeqDataSet(m=48))
design = expand.grid(id=1:3,
phenotype=rep(c("sick","healthy"),each=2),
type = rep(c("A","B"),2))

And you make a group that is a combination of phenotype and type:

design$group = paste(design$phenotype,design$type,sep="_")

Run DESeq2, including id:

dds <- DESeqDataSetFromMatrix(countData=counts, colData=design,
design = ~ id + group)
keep <- rowSums(counts(dds)) >= 10
dds <- dds[keep,]
dds <- DESeq(dds)
res <- results(dds, contrast=c("group", "healthy_A", "sick_A"))

So this compares only healthy_A with sick_A.. hopefully this is what you need... The other option (normal statistical) is to fit a model with interaction, but I think this is easier to explain and interpret.

Another question you might have is why not do two separate analysis for A and B, I think it's sometimes better if you model the variance with all available data, given the limited number of replicates.


How to calculate Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation from protein (amino acid) sequences? I want to code Moreau-Broto autocorrelation, Moran autocorrelation, Geary autocorrelation in Python 3.6 from the scratch for amino acid sequences without using any package. Unfortunately, I have not been able to come across a proper literature which shows how these autocorrelation indices have to be calculated. Can anyone give a clear explanation as to how to calculate it or at least give a link which has got a clear description of how to calculate these?

Rajarshi Guha describes the formulas to calculate Moreau Broto autocorrelation, Moran coefficient and Geary coefficient and has a bibliography.

But why reinvent the wheel? It seems that the Autocorrelation module of the PyBioMed package does already what you want: it is a pure Python implementation without any external dependencies.


Convert FASTA to FASTQ with dummy quality scores I have a FASTA file which I would like to convert into FASTQ format as the tool I want to use my data in requires it in FASTQ format. So dummy quality scores are fine.

Note: I am not using any functionality of the tool that will require the quality score.

A simple Biopython solution:
from Bio import SeqIO

for r in SeqIO.parse("myfile.fa", "fasta"):
    r.letter_annotations["solexa_quality"] = [40] * len(r)
    print(r.format("fastq"), end='')

Example:
$ cat << EOF > myfile.fa
>1
ATG
>2
ATGTAGA
EOF
$ python3 myscript.py
@1
ATG
+
III
@2
ATGTAGA
+
IIIIIII


Subset smaller BAM to contain several thousand rows from multiple chromosomes There are many cases whereby I would like to subset a BAM to create a small file in order to work with (e.g. algorithmic testing, debugging, etc.)

Normally I do the following, which will subset the BAM file.bam and keep the header

samtools view -H file.bam > header.sam
samtools view file.bam | head -n 5000 | cat header.sam - | samtools view -Sb - > file.unique.bam

In this case, I would like 5000 rows in chromosome 1 and 5000 in chromosome 2.

I could first try grepping by individual chromosome, and then combining the two SAMs

e.g. here's complete BAM with grepped chr1 and (incorrect but complete) header

samtools view -H file.bam > header.sam
samtools view file.bam | grep "chr1" | cat header.sam - | samtools view -Sb - > file.unique.bam

but then I have two problems:

(1) I may not be grepping the alignments to chromsome 2---there may exist BAM rows which contain 'chr2' but are not alignments.

(2) I think one must manually edit the header. There's probably no way around this.

Is there an easy way, Bioinformatics SO?

If you're not too hung up on exact numbers like 5000 reads then you can do that with a single samtools command:

samtools view -bo subset.bam -s 123.4 alignments.bam chr1 chr2

That will select 40% (the .4 part) of the reads (123 is a seed, which is convenient for reproducibility). The convenient part of this is that it'll keep mates paired if you have paired-end reads. For 5000 reads per chromosome just change the .4 part to a sufficiently small number.

In general you don't really need to subset the header. Some tools will perform a bit better if you do, but you'll generally get the same results regardless.


What is the aim of insertion codes in the pdb file format? In PDB file we always see insertion codes (ic) in a few proteins. Could anybody please tell me why it has been added in structure and what is the aim of that ???

When I wondered it myself years ago and I found a good explanation in a mailing list archive. I think it was this one, by Frances C. Bernstein:

  The problem with numbering started when people wanted to compare the
  'same' proteins from different species.  They found that there were
  the following possibilities that gave rise to differences:

  More or fewer residues at either end.
  Extra residues at various places within the chain.
  Fewer residues at various places within the chain.
  Different amino acids at the same place.

  Now imagine that residue PHE 195 is very important for the activity of
  the protein in species A.  But in species B it is residue PHE 197 and
  in species C it is PHE 212, because species B and C are not the same
  length as species A.

  Because people felt it was important to preserve the amino acid numbering for 'important' residues and to be able to readily discuss
  and compare the structures from different species, various people
  decided to try to number the proteins from species B, C, etc. to match
  the numbering used for species A.  In doing this, one must have gaps
  (missing numbers) where a sequence is shorter.  But what should one do
  when a sequence is longer?  This is the case where it is necessary to
  insert extra numbers and this is done by using insertion codes.

  Thus the insertion code is an integral part of the residue number and it is improper to ignore that field when using a PDB entry.  You
  must also allow for 'missing' numbers when using a PDB entry.

In the PDB format the insertion code is in column 27, following the sequence number (columns 23-26):

                      nnnnI
ATOM  11918  CZ  PHE D 100      -6.852  76.356 -23.289  1.00107.94           C
ATOM  11919  N   ARG D 100A     -9.676  74.726 -19.958  1.00105.71           N
...
ATOM  11970  CE  MET D 100H     -8.264  83.348 -19.494  1.00107.93           C
ATOM  11971  N   ASP D 101     -11.329  81.237 -14.804  1.00107.41           N

Together, columns 23-27 make a sequence ID, which is a mostly-numeric ID of the
residue position.

Note: It's not a full ID of the residue. When you have a point
mutation (a.k.a. microheterogeneity) you have 2+ residues with partial occupancy at the same position. To fully identify a residue both the sequence ID and residue name are needed.

As of 2020, the PDB (organization) uses PDBx/mmCIF as the primary format. The PDB format is only a legacy one. The presentation of the insertion codes (and many other things) in the mmCIF format are rather confusing, so it's worth explaining.

The mmCIF format introduced double identifiers for chains, residues and atoms: author's ID and "label" ID. Both look the same and are easy to confuse. For example: two chains named B and A by the author are renamed as, respectively, A and B by the wwPDB; both names are used in parallel: one chain is B(A), the other A(B). Poor design choice at every step. The same with the sequence ID. Originally, we had two tags:

_atom_site.auth_seq_id -- author's sequence ID, possibly with an insertion code (100A),
_atom_site.label_seq_id -- numeric-only sequence ID, numbered with respect to the sequence of this molecule (SEQRES).

Then, as I heard, numeric insertion codes were spotted in the wild. (In a database big enough any mistake that is not automatically checked for will happen). This made a problem with translation from the PDB format to mmCIF. The sequence number 100 with IC 1 would be interpreted as sequence number 1001. The numeric ICs must have been remediated in the meantime, because they absent as of 2020. But at that time (10-20 years ago?), instead of correcting the entries, the decision was made to move the insertion code to a separate column (pdbx_PDB_ins_code).
The documentation was never updated, so the auth_seq_id is still documented to store the sequence ID, not just the sequence number.

Here is an example from 3B9F:

                             !!  I                                      nn
ATOM   770  C CD2 . LEU B 2  46  ? -23.416 21.041  -27.401 1.00 19.94 ? 60   LEU H CD2 1
ATOM   771  N N   . TYR B 2  47  A -21.191 16.717  -28.119 1.00 18.73 ? 60   TYR H N   1

The sequence IDs here are 60 and 60A. Since the meaning of the insertion code is not widely known, it landed next to label_seq_id which has nothing to do with it. So again: the insertion code is just a part of the author's sequence ID (originally, in auth_seq_id), and the letter A above (10th field) is simply an extension to the number 60 (17th field).


Comparing gene abundances between metagenomes My workflow until now:
Find fragments of a marker gene in unassembled metagenomes > download and assemble metagenomes > recover the gene neighborhood / gene set of interest
Right now I have a rough estimate of how abundant these genes are by the 'depth' noted on the assembled contigs (I use MEGAHIT for assembly). I was wondering if there's a more thorough/proper way to do this. I would like to compare the abundance of specific genes between a) samples in the same study, and b) different studies. I imagine that the size of the individual metagenomes should be considered in both cases, but point b) might add additional difficulties such as different sequencing techniques. I would appreciate your insights.

I would avoid using assemblies to answer this question, as there's no guarantee that you will be able to assemble your genes of interest; you can however estimate their abundance even if they are relatively rare.
How I understand your question as being one of estimating the abundance of either some specific genes (e.g. butyrate metabolism genes) or all genes in a microbial community across multiple samples for comparative purposes. In other words, not 16S or marker gene analysis for the purposes of estimating organismal abundance, which is a rather different problem (though in that case I would still not use an assembly).
A more standard workflow is:

align metagenomic reads against some existing database of genes annotated appropriately.
estimate the number of reads aligning against some gene or ortholog in some way (using for example KEGG Orthology or similar).
use counts from (2) as input to some statistical procedure, possibly summarizing across functional categories.

Some examples of how this has been done are here, here, here. I am sure that there are more recent/relevant references but I haven't been following the field closely in the last few years.


What's a template switching site? Reading Islam et al. (2011):

  From each transcript, a single read is obtained, corresponding to a template-switching site located preferentially at the 59 end of the mRNA.

By reading this page I think that a template switching site is an additional set of nucleotides added at the 3' end of the newly synthesized cDNA, that serve to anchor a oligo-nucleotide used by the reverse transcriptase to synthesize the second cDNA strand. It is called template switching because reverse transcriptase switches strand upon reaching this set of nucleotides.

Is it correct?

Islam, S. et al. Characterization of the single-cell transcriptional landscape by highly multiplex RNA-seq. genome.cshlp.org
21, 1160–1167 (2011).

From the page you cited:

  During first-strand synthesis, upon reaching the 5’ end of the RNA template, the terminal transferase activity of the MMLV reverse transcriptase adds a few additional nucleotide

That usually are CCC (1st is almost always C and the other two most of the times), see: Zajac et al., 2013

These CCCs are used for priming the Template Switching Oligo (TSO), that has GGG and whatever you want to add to the sequence (for instance UMIs and adapters).

So, answering to your question the polimerase does not switch strands, it switches from mRNA as a template to TSO as template. Because of that is called template switching.

You can see it in the graph from the page you cited. In the template switching step polimerase switch from the mRNA to the TSO and continue with the synthesis of the first cDNA strand using TSO as template.


Build a Pubmed query given long gene list I would like to search the literature using the Pubmed search engine for a set of genes that are associated with a brain region of interest. Unfortunately, there are a lot of them. It's possible there is no solution if there is a character limit on the query. I think it would be easiest to do this with 'nix tools but I'm open to R libraries or Python implementations too.

The list is 908 genes long. It's from the Gensat annotation search engine, with the filters Structure IS olfactory bulb & Expression level IS moderate to strong signal & Expression pattern IS region-specific.

The list looks like, scraped from the dropdown list in the HTML source with a regular expression:

Aatf
Abcb1b
Abhd3
Abhd4
Abhd5
...

The Pubmed query should look like:

(((((((Aatf) OR Abcb1b) OR Abhd3) OR Abhd4) OR Abhd5) OR Abl2) AND olfactory bulb

but with 908 genes separated by OR and parentheses.

Bonus points if you find/write a utility which adds gene synonyms to the query!