STRUCTURE interprets pop ID as locus I'm using structure for the first time and I fail to execute a "job". For some reason, structure interprets the pop ID column as locus.
I have following .txt file

The first column are sample labels, the second column are population IDs and then there are 26 columns containing the locus information (i.e. 13 markers). Thus, in total 28 columns. It is a diploïd organism so each locus consists of 2 consecutive columns (as described in structure manual). In total there are 64 rows (so 64 samples)
I load this .txt file into structure, specifying following parameters:
-number of individuals 64
-Ploidy of data 2
-number of loci 13
-missing data value -9
-special format: data file stores data for individuals in single line
-individual ID for each individual
-putative population origin for each individual
The loaded file looks like:

This looks good. My samples have a separate column, pop ID is a separate column for possible populations an then for each locus there are 2 columns. Not all loci are visible on screenshot, but there are 13.
I want to calculate Delta K and therefore I want to run a job. However when I run it, I get following error:

I get this error for each row/individual. It appears that structure interprets the pop ID values ("Noord","Zuid", "onbekend") as 14th locus, but I have never specified a 14th locus.  I fail to understand why this happens.
When I run the program without pop ID column, it works.
Could anyone please advise me on how to solve this error?

Hmmmm .... you don't actually need pop ID to work for Structure, providing their are no logical errors (you get a wrong answer). The bit that must work is the a priori number of populations in the data set. You are then going to manually (or via a script) feed a range of a priori values into Structure, run each and plot the likelihood. Where the likelihood plateaus thats the number of populations contained within your data. This is the key output of Structure.
Personally, I would perform a Structure run on the test data and remove the pop ID column and check the answers (likelihood -lnL) are identical. Thereafter use an approach such as $F_{ST}$ to interrogate the exact members of each population delineation.

I think your test will work .... The only thing I can think of is the header comprises a space, This format just needs checking against the test data.


How to visualize called narrowPeak files in UCSC Genome browser? I have this file:

and I get this error
I googled and found this solution but it still doesn't work and I get this error

Error File 'GSM2797523_FOXA2_IDR0.02_narrowPeak_try.bed' - line 1 of custom input: variable needs to start with letter '"Peak name=""Somite narrowPeak""" '
Any suggestions?
Many thanks
Found a way!!!
opened it with text editor and not exel as before and added the line:
track type=narrowPeak name="Somite narrowPeak"

Headers should not contain tabs:
track type=narrowPeak name="Somite narrowPeak" description="Somite narrowPeak"

Ensure that you have NO tabs on that line.


Decontaminating RNA seqs for de novo transcriptome assembly and annotation of novel eukaryotes I have raw paired-end RNA-seq reads for two novel eukaryotic species. Some background: the reads represent a copepod (arthropod) species each. The mRNA for each read set was obtained by extracting mRNA from multiple individuals. I am not exactly sure what sequencer was used (some Illumina short read machine as far as I can tell; I am not directly involved in this project).
I'd like to assemble these reads into transcriptomes de novo as no reference genomes are available as far as I know; my hope is to assemble with Trinity and then annotate with Dammit using its default databases (and check for completeness w/ BUSCO).
However, it looks like something went wrong during the sequencing runs (or during sample preparation) because the sequence GC content plots from FastQC look like this:

Here the curve with the tall peak is one sample, and the more flatter curve is the other. There are two closely related curves per sample because FastQC apparently seems to have recognized each read (forward and reverse) as a separate sample.
I first presumed that this is just rRNA contamination (I believe the mRNA was prepared via poly-A selection), and I filtered the reads against all of sortmerna's default databases (which I believe include sequences from SILVA and Rfam among others). Then I passed the sortmerna'd read through FastP to filter out adapters, paired read correction, and low complexity filtering (alongside whatever else FastP does by default). I had ran FastQC before FastP to check if the rRNA filtering did anything and it hadn't, so I also asked FastP to filter by length (--length_required 151) since there were small amounts of reads with lengths shorter than this value which I thought could be the source of the additional peaks in the data.
However, even after sortmerna and FastP together, nothing has really changed:

(There are only two traces/curves here instead of four because the paired reads were apparently recognized properly by FastQC now.)
My plan was to run Khmer to clean the reads a little bit more and then assemble them. I did do a quick assembly with Trinity (with the raw reads) and the BUSCO scores were pretty bad (only ~ 25% of the arthropod BUSCO genes present in the assemblies). I was hoping I'd be able to get Dammit to annotate whatever was left after the cleaning, despite the poor BUSCO prognosis.
Now I'm not entirely sure what to do if I can't clean the reads to begin with, since preventing nonsensical contigs would be absolutely important for a decent annotation.
So this leads me to my questions:

Are these data even salvageable?
Is there anything else I can do to clean the data?
Is this even contamination to begin with that I'm seeing here? (I am presuming some kind of contamination because both samples share that peak at ca. 18% GC content, among other shared features.) Would mapping all the reads against each other and discarding the shared sequences help? If yes, is there a tool I can use for this?
Any other thoughts in this regard?

I searched through the posts here and on Biostars as well as I could but I wasn't really able to find anything helpful in this regard (hence my own post).
I'm very new to all of this, so any inputs, tips, advice, and/or guidance would be most appreciated.

according to your Fastqc files, your reads have a moderate quality. But I am sure there is something to do to improve them.
Try to run fastp with the following options:
fastp -i _1.fq -I _2.fq --detect_adapter_for_pe -l 36 -5 -3 -q 20

and run FastQc again to check to efficacity of the trimming.


PCA vs tSNE in single cell RNA-seq What makes tSNE being the preferred dimensional reduction for visualization in single cell RNA-seq over PCA?

I am aware that tSNE works better at showing local structures and fails to capture global structures of the data.

But I think I don't fully get the reason of why is this an advantage? It offers better resolution? Or better separation of the cells as compared to PCA?

tSNE often offers better visual representation (separation) on such complicated data than PCA. As Micheal pointed out, computing a tSNE embedding over 20.000 gene dimensions is computationally unfeasible, so a number of PCs are normally calculated and these are used as input for calculating the tSNE. They are used in tandem.

As for global vs. local, we are much more interested to see similarity of cells to a limited number of neighbours indicating a celltype and grouping these close together. This is more important than the distance between such celltypes. (assuming the separation in your tSNE is driven by a biologically meaningful factor such as celltype    and not some confounder)

Edit:

Since I just made these images for myself anyway I might as well post them. The first 2 PCs show some separation by celltype, tSNE computed over 100 PCs gets very nice separation.


Need a alternative or more complex version of venn diagram in python for matching dna sequences I am in google colab and I have combined and set up a data frame list of sequences that goes like this
Location ID Sequence
1.1     ........ A ........ AAGAGATA
1.2     ........ A...........   CGGATGAT
1.3     ........ B................   ACGTGATC
1.1     ........ B ........ AAGAGATA
1.2     ........ B...........   CGGATGAT
1.3     ........ C................   ACGTGATC
So I need to make a alternative to Venn diagram because there are about 20+ ID types within a couple hundred rows of list just like this. And some of the sequences repeat or are unique to its id so I need something to be like
CGGATGAT Shows up 4 times in B,C,A, and D
ACGTGATC shows up 1 time in B
AGTGTCTC shows up 5 times in E,A,B,C, and D.
And then have some kind of graphical visualation to show all of these common sequences and where they are located and which are unique. Using python. This might be asking for too much but I'm not sure. I would use something like vennr to make a bigger complex version of this but I don't really know R.

To get around limitations in using Venn diagrams for set overlaps, when there are more than three sets, back around 2013 or so I created something I called an Eulergrid plot (example at the bottom of the page), which an UpsetR plot appears to recapitulate, today.
The Eulergrid code I wrote was in a mix of Perl and R; the UpsetR plot code uses R. There appears to be a Python port of UpsetR that might work, as well.


A Database for Plant Metabolome or Secondary Metabolites? I am searching for a plant metabolome or secondary metabolite database for a college bioinformatics project. Up until now, I have found these two databases, as I am giving their links below;

http://scbt.sastra.edu/pmdb/
https://www.plantcyc.org/https://www.plantcyc.org/

and searched for different scientific articles and reviews;

https://www.frontiersin.org/articles/10.3389/fpls.2013.00073/full
https://www.ncbi.nlm.nih.gov/pubmed/27258266

But only one and updated database which could contain all known plants and all of their known metabolites would be perfect. Could you suggest me any database or any other source with a high coverage for my project?

PhytoMine is an InterMine interface to data from Phytozome, the Plant Comparative Genomics portal of the Department of Energy's Joint Genome Institute, which hosts 81 assembled and annotation genomes, from 78 viridiplantae species.


Why isn't there a standard method to convert GFA to JSON? I'm confused why the GFA format doesn't have a standard JSON parser. It appears to be related to a decision made several years ago, here.

The objection appears to be "JSON isn't necessary, as it's just a simple linked list. Users could write their own". However, I would think this is a useful feature.

There are applications I could think of whereby a JSON would be necessary, or at least very convenient.

(1) Is there a standard way to convert GFA into JSON? Or do users write their own scripts?

(2) Is there a conceptual reason why GFA doesn't convert into JSON?

Short answer: no, there isn't a "standard" method, and there won't be.

Long answer: like TAB-delimited and XML, JSON alone is not a specific format. You have to define a schema to give meanings to data. However, unlike TAB-delimited "format" and XML, there is not an official way to define JSON schema. There have been various attempts but all of them are complicated and none of them are widely adopted. As you can't write down a JSON schema for GFA in the first place, there won't be a "standard" way to convert GFA to JSON. The best you can do is to roll out your own in-house solution, unfortunately.


Extract mapping coverage from GTF files I am doing paired end RNA-seq analysis. I used STAR mapper and then StringTie for getting quantification of my data. After successfully running StringTie, I got .gtf files which will give me information about genes, their coverage along with TPM and FPKM values. It also produced some more .ctab files.

I want to know how what percentage of my reads mapped to reference genome. I know this can be done with BAM or BED files produced after STAR alignment. However, I want to know if this can be done with output files generated by StringTie?

I just want to clarify my question after reading comments and answers I am getting. I want to know if this can be done with .gtf files of StringTie. As I mentioned before, I already know how to achieve this with BAM files. This is little bit weird situation because, I have deleted my BAM files to get more space on my hard drive. So I could just run entire analysis again, but wanted to know if it is possible without redoing whole analysis. Apologies for the confusion.

So consider this, you know how many reads map to each of your genes (1). Does that tell you anything about (2) the number of reads that mapped against non-gene regions (3) the number of reads that did not map against the genome at all?

If you only use the gtf files, you know the number of (1) but not of (2) and (3). To me it sounds like you would want to do something as described in this post.

What it boils down to is:

Finding out how many (unique) reads are mapped, counting pairs once

samtools view -F 0x4 foo.sorted.bam | cut -f 1 | sort | uniq | wc -l

Counting how many total reads you have. Since you have paired-end data the number of forward reads suffice:

cat forward_reads.fastq | grep '@' | wc -l

The percentage is then simply mapped_reads/total_reads


Imputation of dog genotype I am an undergrad working on a set of genotype matrix of dogs, it's in csv file format with information on only the each sample as rows, chromosome and its position as columns, and the rest of the matrix is filled with 0,1,2 and missing values.
So I am wondering what would be the best method to impute the genotype with all the data I get? I looked up Impute2, but it seems to require the reference data(recombination map) which I don't have...
Really appreciate whoever answers my question!!

Impute2 is a great, if dated choice. Another great choice is beagle. The advantage of beagle in this case is that it is not as strict about requiring a recombination map. Without a map, results will be poorer, but perhaps good enough.
Do make sure that your sample size including a reference panel is at least 200 in any case,


Genomics: Downloading genomic data from the NCBI Genbank database I'm taking a genomics course at school, and we have a final project due in a few weeks.
Now, I can write bash scripts to sort and analyze the data, etc. However, I'm not quite sure how to actually get the data from the database. I was hoping to get some help there.

I'd like to get the assemblies of genomes for the species in the Borrelia genus. Is there any way that I can do this using commands like "wget", for example. I'm not even sure what URL I'd have to give the computer to go fetch these. Then once I give it the correct URL, how should I tell it what file in that directory to download?

I'd greatly appreciate some help on this issue.

For downloading data from commandline use the E-utilities of NCBI.


How to use busco installed through conda? I have installed busco on my Mac computer using conda install --channel bioconda busco=3.0.2 into an environment called busco3_0_2. However, when I type busco, the command is not found.

Anyone knows what the name of the binary is? Usually when I install something with conda, the name of the command is the same as the name of the package. In this case, it does not seem to be so.

Note that the package installation finished successfully.

From the README, it looks like the command-line entry point is not busco but run_BUSCO.py. This seems to agree with the bioconda recipe, where run_BUSCO.py and run_busco are listed among the commands.


What command/invocation is used to generate NCBI 16SMicrobial blastdb I'm looking for the exact invocation used to generate the 16SMicrobial database that you can download from here:

https://ftp.ncbi.nlm.nih.gov/blast/db/

I'm hoping to create the same type of blastdb with the same type of metadata with custom sequences.

Platform isn't an issue, but let's say on ubuntu 14.04 or 16.04.

I would like to replicate the creation of the database as closely as possible.
 The most important feature is the taxonomic information as can be seen here:

  The databases on the FTP site contain taxonomic information for each
  sequence, include the identifier indices for lookups, and can be up to
  four times smaller than the FASTA. The original FASTA can be generated
  from the BLAST database using blastdbcmd

Creation of a blastdb using makeblastdb from a set of a fasta sequences is not an issue and can be achieved via:

makeblastdb -in <your_file.fasta> -dbtype nucl -out <database_name>

My question is specifically about the invocation NCBI uses to add the metadata that is present in the NCBI's 16SMicrobial blast database as I am keen to make sure I have replicated the process as closely as possible.

Running blastdbcmd -db foo -info provides a little information but I haven't seen anything which will report exactly how a blastdb was created.
A blastdb readme suggests that only the parse_seqids option has been added to the standard parameters.

For those from NCBI, the following makeblastdb commands are recommended:

For nucleotide fasta file:   makeblastdb -in input_db -dbtype nucl -parse_seqids
For protein fasta file:      makeblastdb -in input_db -dbtype prot -parse_seqids

The -parse_seqids requires NCBI fasta headers in the correct format (spec can be found here).


Repeatscout installation problem I've installed Repeatscout from https://bix.ucsd.edu/repeatscout/.
I've also installed TRF and nseg required by the program.
I followed the steps:

build_lmer_table -l 16 -sequence ~/file_fasta/od_1000.fa -freq od_1000.freq

 RepeatScout -sequence od_1000.fa -output od_1000_repeats.fa -freq od_1000.freq -l 16

and the I put this command:

cat od_1000_repeats.fa| filter-stage-1.prl >od_1000_repeats_filt1.fa

but this message appears:

bash: filter-stage-1.prl: command not found...

Why? It is all correctly installed, I can't understand the problem

If you're in the RepeatScout directory you need to execute the Perl script with a ./ prefix, like so.

cat od_1000_repeats.fa | ./filter-stage-1.prl > od_1000_repeats_filt1.fa

If you're not in the RepeatScout directory, you need to provide the full relative or absolute path when executing the script. Alternatively, you can put the RepeatScout directory in your $PATH environmental variable.

These are all standard UNIX/Linux concepts that many bioinformaticians take for granted, for better and for worse. If you're lost in all the terminology, a crash course in UNIX basics is probably where you would want to start.


Generic HMM solvers in bioinformatics? Hidden Markov models (HMMs) are used extensively in bioinformatics, and have been adapted for gene prediction, protein family classification, and a variety of other problems. Indeed, the treatise by Durbin, Eddy and colleagues is one of the defining volumes in this field.

Although the details of each of these different applications of HMMs differ, the core mathematical model remains unchanged, and there are efficient algorithms for computing the probability of the observed sequence given the model, or (perhaps more useful) the most likely hidden sequence given the sequence of observed states.

Accordingly, it seems plausible that there could be a generic software library for solving HMMs. As far as I can tell that's not the case, and most bioinformaticians end up writing HMMs from scratch. Perhaps there's a good reason for this? (Aside from the obvious fact that it's already difficult, nigh impossible, to get funding to build and provide long-term support for open source science software. Academic pressures incentivize building a new tool that you can publish a paper on much more than building on and extending existing tools.)

Do any generic HMM solver libraries exist? If so, would this be tempting enough for bioinformaticians to use rather than writing their own from scratch?

I would also recommend to take a look at pomegranate, a nice Python package for probabilistic graphical models. It includes solvers for HMMs and much more. Under the hood it uses Cythonised code, so it's also quite fast.


How do I estimate branch length of a maximum parsimony phylogeny in R? I am trying to estimate a phylogenetic tree in R with ape and phangorn. With pratchet() and optim.parsimony(), I can estimate a phylogenetic tree, but branch lengths are not estimated. Is there a way to estimate a phylogeny with branch lengths with the following made-up data:

x <- data.frame("L1" = c(0,0,1,0,0,0,1,1,1,1,0,0,0,0,1), "L2" = c(0,1,1,0,1,1,0,1,1,0,1,1,0,1,1),"L3" = c(1,1,1,0,0,0,1,1,1,0,0,0,1,1,1), "L4"=c(0,0,1,1,0,0,0,0,1,1,0,0,0,0,1))

@Nathan S. Watson-Haigh's answer does give branch lengths, but NOT parsimony branch lengths, which may be what the OP wanted.  nnls.phylo is giving distances calculated from a pairwise distance matrix. To get the parsimony branch lengths (i.e. number of changes mapped on to each branch), score the tree against the dataset using acctran:

treeRatchet <- pratchet(primates)
treeRatchet.BL <- acctran(treeRatchet,primates)
plot(treeRatchet.BL,main='Phylogeny (Parsimony branch lengths)')

Note that the branch lengths in the parsimony phylogram are 'chunkier' than those in the nnls.phylo phylogram -- that's because

they are all integers, and
they aren't corrected with reference to a model of substitution that considers multiple hits or base composition.


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

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

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

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

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

And I am getting an error:

  Error in .checkedCall(cxx_colsum_subset, exprs_mat, is_feature_control) :
    subset indices out of range

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

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


Mapping UniProt id to Entrez id How can we get mapping from UniProt id to Entrez ID? I have UniProt id, but want to convert those ids to Entrez or Ensembl id.

The Retrieve/ID mapping tool from Uniprot allows you to convert from Uniprot to entrez ID and vice-versa. It is available here. Just select the appropriate values in the 2.Selection options panel.


Stable download URLs One big problem that I'm regularly facing is that URLs for downloading Bioinformatics data (e.g., RefSeq releases or NCBI genome releases) disappear.

Does anyone have any good solution for this?

The Persistent uniform resource locator or PURL is one such solution, these are designed to be a bit more robust than permalinks in so much as they are supposed to survive the change of domain name.  The bio ontology community already use them http://purl.bioontology.org/docs/index.html


Why do NEBNext indexing primers have sequence between the p5 oligo and index? In a previous post I asked Why do NEB adapters have non-complementary sequence?

Since then, I realized that there is some other sequence in the p5 indexing primer, as well as in the p7 indexing primer.

Here is a diagram of the NEBNext protocol. The parts that I am confused about are in step 5.

My first question is about the p5 index. Why is there a sequence 5'GATCTACAC 3' between the p5 flowcell oligo and the index?

My second question is about both indices. Why is there no gap between the p7 index and the flowcell annealing oligo?

For reference, I fetched these indexing primer sequences from NEB on page 18 and 19 of NEBNext Multiplex Oligos for Illumina instruction manual.

Question 1: The additional sequences is needed because that is complimentary to the P5 sequence anchored to the flow cell. It is also the site for priming the Index 2 read on a MiSeq. The additional sequence is also needed for the read 1 primer in the cartridge to anneal to the correct place for the molecules.

Question 2: The "GT" sequence you're referring to is necessary for the Index 1 and Read 2 primer to anneal to the correct place.


Is there any way of using biopython to write Swissprot files? I have the following mwe for filtering a Swissprot file based on a certain feature, in this case, transmembrane regions.

from Bio import SeqIO
records=[]
for record in SeqIO.parse("Input.txt", "swiss"):
    transmembrane_protein=False
    print record.id
    for i, feature in enumerate(record.features):
        if feature.type == "TRANSMEM":
            transmembrane_protein=True
    if transmembrane_protein==True:
        records.append(record)
SeqIO.write(records, "Output.txt", "swiss")

The script works when SeqIO.write(records, "Output.txt", "swiss") becomes SeqIO.write(records, "Output.txt", "fasta")

However, such a method is not yet supported.

ValueError: Reading format 'swiss' is supported, but not writing

From the docs, I see that writing is not supported for swiss:

  Note that while Bio.SeqIO can read all the above file formats, it cannot write to all of them.

Is there any unofficial way of using biopython/python to write swissprot files from parsed swissprot files?

Using SeqIO.index rather than SeqIO.parse lets you read all the records into a dict, from which you can then extract the raw entry:

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
for key in record_dict:
    print(record_dict.get_raw(key).decode())

Now you should be able to apply your test for a transmembrane protein to each entry, and write out only the ones you want to keep.

from Bio import SeqIO

record_dict = SeqIO.index('Input.txt', 'swiss')
out = open('Output.txt', 'wb')
for key in record_dict:
    record = record_dict[key]
    for i, feature in enumerate(record.features):
    if feature.type == "TRANSMEM":
        out.write(record_dict.get_raw(key).decode())


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

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


Why PHMM is used For Protein Analysis I am interested on how profile hidden markov model is used for protein analysis in bioinformatics. Can you guys explain why phmm is more preferable method comparing to simple MSA(Multiple Sequence Alignment) scoring or simple MSA probability scoring? For example we have three MSA sequences as based model and a target to be scored.

MSA Sequences

A B C B C A
A B B C C A
A B C B A A

Target Sequence

A B B B C A

By MSA scoring we can count score for target sequence comparing to a profile MSA. For example above target sequence protein column will counted as match(A), match(B), mismatch(B-majorityMSA(C)), match(B-majorityMSA(B)), match(C-majorityMSA(C)), match(A). Just give score for match and mismatch and we can get total score for a target sequence comparing to a profile MSA sequences.

By MSA probability scoring we can count probability of each protein in each column so for example above
P(target_seq similar to MSA profile) = P(A in MSA sequences)*P(B in MSA sequences)*.... *P(A in MSA sequences) = 1*1*(1/3)*(2/3)*(2/3)*1 = 4/27

Isn't MSA scoring like above examples powerful enough?

Thanks a lot !!

The simple probabilistic scoring scheme you describe to query a target sequence against a (protein) family represented as a multiple sequence alignment (MSA) is actually similar in essence to a Position Weight Matrix (PWM), although PWMs are more refined. PWMs are powerful rapid tools often used for motif analyses.

However, your approach (and PWMs) don't handle insertions or deletions (indels), and you can't easily add these into your model (for example, consider that 2 1bp indels should not score the same as 1 2bp indel, and indels could occur within your MSA as well as in your target sequence). You could try and add some arbitrary gap opening and extension penalties somehow (like BLAST) but then you'd abandon your probabilistic model! By contrast, pHMMs have a fully fledged probabilistic indel model with insertion and deletion states.

Finally, of practical note, there are some excellent implementations of pHMMs for protein sequence analysis, such as HMMER3 and HHSuite.


How to merge back the processes by sample after splitting by chromosome in Nextflow? I have a pipeline that has an input list of samples. I would like to split those samples by chromosomes, process each chromosome apart, and then merge them back by sample name.
Here's an example:
samplesList = Channel.from(["sample1", "sample2"])

// A process that splits the samples by chromosome as an example
// In my actual case, I have a VCF that is split by `snpSift Split` command
// and generates files like [sample]_[chr].vcf
// In the example below, I only gave sample1
// But take into consideration that I have multiple samples

process bar {
   echo true

  output:
    file("*.tsv") into bar_output

  script:
  """
  echo 1 > sample1_chr1.tsv
  echo 2 > sample1_chr2.tsv
  echo 3 > sample1_chr3.tsv
  echo 4 > sample1_chr4.tsv
  echo 5 > sample1_chr5.tsv
  """
}

// Process each chromosome apart for each sample
// In my case, I just work some magic
// On the file and then output another file
// This means the output of the process foo
// is only ONE chromosome per sample
process foo {
  input:
    each x from bar_output  // to treat each file separately from the channel
  output:
    file("*.tsv") into xasxa
  script:
    """
    echo \\\$(cat $x)\\\$(cat $x) > sample1_chr\\\$(cat $x)_2.tsv
    """
}

// What is the next process to group back them together by sample name?
// I do not wish to do a collect() on channel foo because it will force
// all the processes of foo to finish first.
// In this example, foo would emit 5 different channels of chromosome
// for sample1 (imagine also there are a bunch of other samples included)
// This process should merge sample1 only when all chromosomes are processed
// independently of other samples.
process mergeChromosomesBySamples {
   input:
     ???
}

As the example above, the process mergeChromosomesBySamples should return something similar/close to
["sample1", ["/path/to/workdir/sample1_chr1.tsv",
             "/path/to/workdir/sample1_chr2.tsv",
             "/path/to/workdir/sample1_chr3.tsv",
             "/path/to/workdir/sample1_chr4.tsv",
             "/path/to/workdir/sample1_chr5.tsv"]
]

EDIT1: more explication for the process mergeChromosomesBySamples given in the comments and as follows ...
The channel xasxa processes each chromosome apart, this means that the output of the channel xasxa contains only one chromosome of the sample like "/path/to/workdir/sample1_chr1.tsv" only. My goal is then to merge the processes out xasxa for each sample together. Sorry if that wasn't clear in my explanation.
With regards .collect() on xasxa I do not wish to use collect() in this case because it would have to force all the processes of foo to finish before. I want them to work independently for each sample.
Where each sample has his own VCF generated and I process them independantly.
EDIT2: Here's a graph of what I am trying to achieve:

If I use .collect(), then all the processing on the chromosomes of all the samples have to finish first before being able to merge them, which is something I do not want.
Thank you for your help.

Use the groupBy operator. If the content of channel xasxa is:
["/path/to/workdir/sample1_chr1.tsv", "/path/to/workdir/sample1_chr2.tsv", "/path/to/workdir/sample2_chr1.tsv", "/path/to/workdir/sample2_chr2.tsv
You could probably do:
xasxa.groupBy{ it -> it.fileName.toString().split("_")[0] }
The result would be a channel with:
[sample1:[/path/to/workdir/sample1_chr1.tsv, /path/to/workdir/sample1_chr2.tsv, sample2:[/path/to/workdir/sample2_chr1.tsv, /path/to/workdir/sample2_chr1.tsv]]
This breaks if you have an underscore in your samplename before _chr*. More robust would be to initialise all channels with a sample value, so that each channel is in the form of [samplename, chromosome]. Then all you need is a groupTuple on your channel.


Find intervals that genes fall within their range I am trying to run some frequency based stats to identify selective sweep in my study system using Rpackage "PopGenome".

To proceed with it I have split my genome data (whole genome sequencing) into 5MB chunks

cregions <- character(31)
cregions[1] <- "1-5000000"
cregions[2] <- "5000001-10000000"
cregions[3] <- "10000001-15000000"
cregions[4] <- "15000001-20000000"
cregions[5] <- "20000001-25000000"
cregions[6] <- "25000001-30000000"
cregions[7] <- "30000001-35000000"
....

I have a bunch of genes and I want to know each of my genes falls within the range of which intervals in the above list.

these are my genes:

chrom start stop gene_id
Chr01 107310787 107323633 gene2748
Chr01 113789595 113794851 gene2894
Chr01 159041389 159068535 gene4806
Chr02 2742629 2758060 gene4846

I'd like to get an out put like this

gene2748  window_2
gene2894  window_3
...

Any idea how I can do it in R?

The GenomicRanges package is a great place to start, something like this:

library(GenomicRanges)

# make a GenomicRanges set corresponding to the genes
genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom")

# turn your windows into a data frame and make a GenomicRanges set from that
windows = ...
windows.gr = makeGRangesFromDataFrame(df = windows ...

# find overlaps between the two sets
olaps = findOverlaps(genes.gr, windows.gr)

# output relevant fields
cat(paste(genes[queryHits(olaps),"gene_id"], windows[subjectHits(olaps),"window_id"]))


CIBERSORT runtime error eval failed I was running CIBERSORT but caught this error right after it started permutation:

org.rosuda.REngine.Rserve.RserveException: eval failed, request status: error code: 127
88

Ironically, I found that this problem was caused by the fact that the tables are too large!
They really should have designed the correct error messages!


install bowtie2 from sources cannot find -ltbb I am trying to install bowtie2 aligner from sources without the root access.
Bowtie2 needs tbb package to be installed and it is recommended to install oneTBB.
# install one TBB
mkdir -p SHOME/soft
cd $HOME/soft
git clone https://github.com/oneapi-src/oneTBB.git
cd oneTBB
make

Then I add C_INCLUDE_PATH to .bashrc:
echo "export C_INCLUDE_PATH=$HOME/soft/oneTBB/include" >> ~./bashrc

and continue to install bowtie2:
# install bowtie2
mkdir -p $HOME/soft/bowtie/source
cd $HOME/soft/bowtie/source
wget https://sourceforge.net/projects/bowtie-bio/files/bowtie2/2.4.1/bowtie2-2.4.1-source.zip
unzip bowtie2-2.4.1-source.zip -d ..
cd ../bowtie2-2.4.1
make

TBB library is recognized but in the end I get this error:
/usr/bin/ld: cannot find -ltbb
/usr/bin/ld: cannot find -ltbbmalloc_proxy
collect2: error: ld returned 1 exit status

ltbb and ltbballoc_proxy are not in TBB but libtbb and libtbballoc_proxy instead.
I also tryed to add variables to .bashrc:
TBB_LIB_PATH=$HOME/soft/oneTBB
TBB_INCLUDE_PATH=$HOME/soft/oneTBB/include
TBB_RELEASE_DIR=$HOME/soft/oneTBB/build/home/suvar/apps/oneTBB/build/linux_intel64_gcc_cc9.3.0_libc2.31_kernel5.4.0_release

but seems it does not help.
I need help to complete the installation of bowtie2 on our server.

You need to give the path to the source files for the compiler and the full path to the compiled libraries (all parent directories included) for the linker.  You also need to use either CPATH or CPLUS_INCLUDE_PATH instead of C_INCLUDE_PATH to point to the oneTBB source, since bowtie2 is C++ and C_INCLUDE_PATH will only have an effect for C (see Environment Variables Affecting GCC in the GCC documentation for more info.)
For example like this, inside bowtie2-2.4.1, without making any lasting changes to the shell's setup in this case:
$ ONETBB=$HOME/soft/oneTBB
$ export LIBRARY_PATH=$ONETBB/build/linux_intel64_gcc_cc5.4.0_libc2.23_kernel4.4.0_release
$ export LD_LIBRARY_PATH=$LIBRARY_PATH
$ export CPATH=$ONETBB/include
$ make

After it's built you'll still need the library path to oneTBB defined (via LD_LIBRARY_PATH) so bowtie2 can find it at run-time.
But if you don't have a specific reason to compile from source I'd recommend avoiding it.  You could go with what h.mon suggested and use an existing binary, or use a package manger to install it like zorbax suggested.  Conda has a package for bowtie2 in the bioconda channel, up through and beyond version 2.4.1, and you can use that as a package manager without root access.


server for finding kmers in set of sequences Is there a server/website somewhere where I can submit a list of DNA/RNA sequences and find the list of kmers hits and organisms where it's found? I checked the kraken website, but they don't have a webserver for it.

I am specifically looking for kmers, as a complement of just doing BLAST on the NCBI nt database or similar.

EDIT: So far I have seen:

the Kraken app on Illumina BaseSpace, but it's not trivial to upload FASTQ data that hasn't been produced from an Illumina instrument.
the CLARK trio of apps in InsideDNA, tried once with default parameters but produced empty results.

One Codex is powerful and user friendly, and works in a similar way to Kraken underneath
https://www.onecodex.com/

The Kaiju web server might be of interest and is more sensitive than exclusively k-mer based methods
http://kaiju.binf.ku.dk/server

The title of this question is misleading – it sounds like you're actually looking for a webserver for taxonomic classification using k-mers.


Microarray Meta-Analysis and Cross-Platform I have a microarray dataset with Agilent platform, all of the samples are cancerous, I need normal or control samples to compare with cancerous samples, can I combine control samples from another platform (for example Illumina) to Agilent platform? and then remove the batch effect (by ComBat function in the sva package).
I would be grateful if you could guide me to do it.

In short the answer is no because you have to standardise between the arrays. Thus both arrays would need to run exactly the same control sample and you would then perform statistical standardisation (its a formal term), which is a varience measure where you look at the variation around the 'mean' for all samples in the array and here the 'mean' is the standard sample. If that can be done for both arrays at that point they can be combined.
The problem is intensity per run varies and you have to control against this

In this scenario you  can certainly combine samples by assessing whether the control samples between runs approximates to identity. What you could use is a uniform distribution to assess this, variance must MUCH smaller than the mean for the STANDARD sample(s). I don't really do uniform distribution statistics, so I'd need to research this, because 'uniform' means identity to standard deviations against the uniform distribution is a bit meaningless. However again the uniform distribution is denoted by a very variance in comparison to the mean. The statistic of dividing mean by standard deviation/variance from memory is called the Index of Association (I could be wrong). If this value is at an extreme you're good to go <<<1, if it is remotely close to 1 you're in trouble, if it is >1 you need to use standardisation.
Presumably you have several standards to ensure a mean is meaningful. A T-test and non-parametric equivalents is not stringent enough. If there is a difference assess a correction factor and rerun the uniform distribution check. The same correction will be applied to all samples.
The undisputed method to combine is standardisation, but accept in this scenario the differences between runs is potentially very small, so here what you are attempting is essentially to combine raw data or at least keep one of the runs as pure raw data.


visualisation of genome alignment I was asked to check the synteny of some genes in two genome assemblies of the same species (PacBio & Illumina). I was given two scaffolds couples (so 4 scaffolds total) in the Illumina genome and asked if they were assembled together in the PacBio assembly.

So, I proceeded with a DNA-DNA alignment with minimap2 of the two Illumina scaffolds on the PacBio whole genome. But I am not sure on how to directly use the sam file for the question I was asked. So I tried to visualize this alignment using Mauve. It worked for the first couple, but it failed (blank page) for the second couple of scaffolds.

Surely, Mauve is not the best option out, so can you suggest an alternative? This is a chordata genome.

After a bit of digging, finally found exactly what I was looking for. I'll share here in case someone else was looking for the same thing.

I used a standard Nucleotide BLAST (blastn), coupled with MUMmer. Short tutorial.

Assuming that xn & yn are the couple of scaffolds I want to check if they are assembled together in the Pacbio assembly (a multi-fasta file):

1) first, I blasted the pair of Illumina scaffolds (as separate queries) to the PacBio genome (as subject):

blastn -query ../illumina_asm/scaffold_x1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_x1; blastn -query ../illumina_asm/scaffold_y1.fa -subject ../pacbio_asm/genome.fa > out_scaffold_y1

The candidates PacBio scaffolds shall appear in the top blastn hits in both out_scaffold_x1 & out_scaffold_y1.

Let's assume we found one common scaffold only - corresponding to the fasta entry pacbio_scaffold in the PacBio genome.

2) extract the sequence of the target scaffold from the PacBio assembly, using samtools faidx:

samtools faidx ../pacbio_asm/genome.fa; samtools faidx ../pacbio_asm/genome.fa 'pacbio_scaffold' > pacbio_scaffold.fa

3) concatenate the Illumina scaffolds to analyze:

cat ../illumina_asm/scaffold_x1.fa ../illumina_asm/scaffold_y1.fa > illumina_scaffolds.fa

4) compare the Illumina scaffolds to the PacBio assembly using nucmer from MUMmer:

nucmer --prefix=pacbio_scaffold pacbio_scaffold.fa illumina_scaffolds.fa

this will generate a pacbio_scaffold.delta file, used in the next step

5) filter the delta file to keep only the most meaningful alignments file using delta-filter from MUMmer (avoiding confusing messy plots later):

delta-filter -m pacbio_scaffold.delta > pacbio_scaffold.delta.m

6) use mummerplot from MUMmer to achieve a graph:

mummerplot --png pacbio_scaffold.delta.m -R pacbio_scaffold.fa -Q illumina_scaffolds.fa --prefix=pacbio_scaffold -large -layout

On the x-axis, the PacBio scaffold they have in common. On the y-axis, the two Illumina scaffolds (separated by the horizontal grey abline). Colors correspond to the orientation.

Example of Illumina scaffolds linked in the PacBio assembly:


Scatter tasks around predefined batches inside a WDL workflow Suppose I have two arrays Array[File] X and Array[String] Y both of length N, with Y being a sort of R-like factor that defines batch membership for X, for a total of M batches with M < N. A valid example would be:

Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

In which case N=5 and M=2. I want to scatter M tasks so that each task takes only the files corresponding to the same batch. In the previous example, one task would take as input:

Array[File] X =["file1", "file2", "file4"]
String Y = "A"

And the other task would take as input:

Array[File] X =["file3", "file5"]
String Y = "B"

Is there an easy way to do this within a WDL workflow in the general case? I understand that there are less elegant ways to do this as its own task but, as this seems like it could be a general purpose scatter/gather problem, an easy solution might be of general interest to the community.

Here's an alternative to using a custom task, albeit still rather unappetizing:

workflow groupby {
  input {
    File? nil
  }

  Array[File] X = ["file1", "file2", "file3", "file4", "file5"]
  Array[String] Y  =  ["A",     "A",     "B",     "A",     "B"]

  scatter (y in ["A", "B"]) {
    scatter (p in zip(X,Y)) {
      File? x_or_nil = if (p.right == y) then p.left else nil
    }
    Pair[String,Array[File]] group = (y, select_all(x_or_nil))
  }

  output {
    Array[Pair[String,Array[File]]] ans = group
  }
}

(It'll be slightly cleaner with the upcoming None literal instead of dummy nil)

FWIW, here's a lengthy discussion about the idea of allowing developers to define custom library functions which just run a small task that manipulates data structures (hinting that it needn't necessarily offload to a worker node, with the huge latency that typically entails).


What is the right way of calculating a Phred score by hand? i am trying to calculate mean Phred scores for my sequencing data, but i feel not very comfortable about it. There are actually two ways of calculating. (I just use an existing sample)
giving: 3 reads with a mean phred score of 20, 10 and 3 (meaning Error rates of 1%, 10% and 50%, respectively)
Way 1 to calculate: Just take the arithmetic mean of the q scores:
(20+10+3)/3 = 11
Way 2:  Take the arithmetic mean of the error rates:
(0,5+0,1+0,01)/3 ~= 0.2
And then calculate the q score by
-10*log10(0.2) = 6.99
(https://gigabaseorgigabyte.wordpress.com/2017/06/26/averaging-basecall-quality-scores-the-right-way/)
The Author of this Blog claims, that Way 2 is the correct version, but i still do not get why.
Is somebody able to answer this question?
Thanks a lot

Not sure if my explanation is any good but let's try...
First, let's convince ourselves that converting phred to probability is the right thing to do (i.e. your Way 2).
Imagine a read of length 10000 bp with half of the bases with Q=10 (or P=0.1) and the other half with Q=20 (or P=0.01). The number of wrong calls you expect is therefore 550:
n = 10000
(n * 0.01 * 0.5) + (n * 0.1 * 0.5) # -> 550 wrong calls

So the overall error rate is 550/10000 = 0.055 which corresponds to Q=12.6. If you follow your Way 2 you do get this answer:
((n * 0.5 * 0.01) + (n * 0.5 * 0.1))/n = 0.055
-10*log10(0.055) = 12.6

If you follow your Way 1 you get Q=15, which means an overoptimistic error rate of ~0.03 or ~316 wrong calls. Why is that?
Because phred is on a non-linear, logarithmic, scale: a unit change  towards the high values is less important than a unit change towards the lower end. For example, across 10000 bases a change from Q=30 to Q=31 means a change from 10 to ~8 wrong calls - a difference of 2 only. But a change from  Q=10 to Q=11 means 1000 to ~800 wrong calls: the difference is a substantial ~200 calls.
Having said this, the straight average of phred scores is not "wrong". It measures the overall phred score but it does not reflect the expected number of wrong calls, which is probably what you want in this case.


Map gene IDs to protein structure if exist I want to map some gene ids to their proteins structures if they exist.
The genes are rows of this CSV. There are other information that may help also such as Product description but that may lead to anbiguities. I am looking for a programmatic way, biopython or pdb API?

Not sure if you want the structure exactly or a handle of some sort, but the best way i know to get ALL your 5000 identifiers to map to uniprot/pdb is uniprot's mapping service.
You can either copy your ids straight from excel or parse with python ex.:

import pandas as pd
import os

file_name = 'aap7847_Zhang_SM_Table-S5.xlsx'
sheetname = 'Table S5'
df = pd.read_excel(file_name, sheet_name=sheetname, skiprows=1)
ids = list(df['Gene_ID'])

name_of_file_to_save_to = "geneIds.txt"
with open(name_of_file_to_save_to, "w") as newf:
    for id in ids:
        print("{}\n".format(id), file=newf)
    newf.close()

You can then save the Uniprot's response and either add it as a column to your file or map to other parts of Elixir ecosystem with SIFTS, for example, as @Matteo Ferla says.

In the case that you maybe need a single id's information, you could hit their API too, although this wouldn't work for bulk or at least didn't work for me. You can map the uniprot accession you get from this to PDB or whatever else you need!

import requests
def get_uniprot_accession(geneid, format='list'):
    resp = requests.get(
        "https://www.uniprot.org/uniprot/?query={}&sort=score&format={}".format(geneid, format)).text
    if format == 'list':
        resp = resp.strip()
    return resp


bed file with N regions of GRCh38 reference? How do I obtain a bed file with the list of N-regions of GRCh38 reference? This is, the regions where the sequence is a stretch of Ns.

# if you have seqtk installed, skip the following two lines
git clone https://github.com/lh3/seqtk
cd seqtk && make
# the main command line
./seqtk cutN -gp10000000 -n1 hg38.fa > hg38-N.bed

Option -n sets the min stretch length. Just use -p as is. It is a bit complicated to explain what it is doing.


Why most aligners do not output the "X" CIGAR operation? As I read the SAM spec, the "X" CIGAR operator represents a mismatch. This seems useful as we can know where are the mismatches without looking at the reference genome. However, many popular aligners such as BWA do not output "X". Why do they omit "X"?

The SAM format originally had only M, I, D, N, S, H, and P CIGAR operators.  See the original SAM specification (if you can view Apple Pages documents) and Table 1 in The Sequence Alignment/Map format and SAMtools (Li et al, 2009).  This was in line with previous tools using CIGAR strings, notably exonerate which introduced them with just the M, I, and D operators.

BWA-backtrack was written contemporaneously with the SAM format in 2008 and 2009 (and published in May 2009).  Its ChangeLog shows that it was outputting mismatch information in an MD tag from January 2009 and that MD was defined in the SAM specification of the time (and that the tag value's syntax was somewhat in flux in February 2009).  The MD tag is also described in that early v1.0 Pages-formatted SAM specification.

The = and X CIGAR operators were introduced later in SAM v1.3 as a result of this lengthy samtools-devel mailing list thread.  The characters used for the operators and the initial implementations were essentially in place by November 2009.

Since then, the =/X operators have not really taken over from M, as you've seen.  There's probably a number of factors contributing to this:

Their later introduction, long after the basic M/I/D operators and MD tag were well-established;
Conceivably their being specific to SAM and unavailable in other CIGAR flavours;
The MD tag still provides more information — X doesn't tell you what the mismatched reference bases were.


Interpreting GWAS results with different settings I did a bunch of GWAS analysis (linear model without covariates) with applying different quality controls. How to choose the optimal settings when filtering for minor allele frequency (maf), Hardy-Weinberg equilibrium (HW) and missing call rates (mcr)? Are there any rules? Here are some examples:
Mahnhatten and QQ plot with settings maf=0.02, HW=1e-10, mcr = 0.2:

Manhattan and QQ plot with settings maf=0.01, HW=0.0001, mcr = 0.01:

Mahnhatten and QQ plot with settings maf=0.02, HW=1e-10, mcr = 0.05:

Based on the QQ plots, filter out variants with p > 0.01, then display the beta statistic (or whatever other statistic you're using to demonstrate association). If you want to incorporate the p-value in the displayed results, show the minimum test statistic from a 95% confidence interval. Here's an attempt I made at doing that with some GWAS summary statistics relating to household income (something that wouldn't be expected to have a heritable component):
https://twitter.com/gringene_bio/status/1207617723586371584
After looking at the test statistic, you may be able to make better decisions about thresholding with the other values. One thing you could try is to incorporate those values into a scatter plot that shows more than two dimensions. The other statistics could be represented as size, colour, or Z location on a 3D scatter plot.
P-values should not be overemphasised, used for ranking, or displayed on their own. See the American Statistical Association's report on the p-value for more information:
https://amstat.tandfonline.com/doi/full/10.1080/00031305.2016.1154108#_i30
As far as I'm aware, the p-values as calculated by GWAS software represent uncertainty in the test statistic (i.e. t, or beta), and not the probability of association. If you want to calculate that association probability, carry out bootstrap sub-sampling of the individuals and then determine the proportion of sub-samples in which the tested variant appears in, say, the highest 5% of variants, ranked by the test statistic. Instead of filtering on HWE, MAF and MCR, you could try filtering on the bootstrap proportion instead, which should incorporate all statistics that introduce uncertainty into the dataset.
More details and figures can be found on this poster:
https://f1000research.com/posters/5-2190
If you're going to be generating a PCA on these data, please consider labelling the cases and controls (or low/high for quantitative traits) on the PCA to verify that the groups look similar in terms of overall population structure. Hopefully you won't find any blobs in the PCA that have only one of the two groups in them, but if you do, they should be dealt with first. Any differences in population structure (even very small, subtle differences) will dominate the GWAS results because structure-associated signal will cover a huge region of the genome in comparison to effects from a few (or a few hundred) SNPs. Here are some papers I've found insightful regarding this issue:

https://doi.org/10.1006/tpbi.2001.1543
https://doi.org/10.1038/s41467-018-05257-7


Downloading URL using AWK for fixed fields I have a text file as follows

a   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz
b   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz
c   k   testis  adult   male    8 week  rRNA    https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz

The file is tab delimited and has 8 columns with 8th column being the url link  I need to download the files one by one using wget and changing the name of the downloaded file with -O option in wget. The argument I want to give to the -O option is combination col1 and col2 separated by _ so that for first line I get the downloaded file as a_k.fastq.gz and so on. For example, the first line of the file above should give:

wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz

Following the post from here I realized looping through the file is slow and sometimes difficult. How can I download these files using wget with awk or cut as the tool for fixed fields

awk '{print "wget " $9 " -O " $1 "_" $2 ".fastq.gz"}' tmp

will give you

wget https://www.encodeproject.org/files/ENCFF001JPX/@@download/ENCFF001JPX.fastq.gz -O a_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JPZ/@@download/ENCFF001JPZ.fastq.gz -O b_k.fastq.gz
wget https://www.encodeproject.org/files/ENCFF001JQB/@@download/ENCFF001JQB.fastq.gz -O c_k.fastq.gz

If you add a | bash to the command it will then execute all the wget commands.


Automating sequence download from Wormbase I'm working on the C. elegans model organism, for which there exists an online biological information database: Wormbase.

I'm able to manually search for a transcript (say "K06C4.12") and this gets me the following URL:
http://www.wormbase.org/species/c_elegans/cds/K06C4.12

Say I have this transcript ID: "F08G2.2". I guess that I can form a valid similar URL, as follows:
http://www.wormbase.org/species/c_elegans/cds/F08G2.2

Indeed, it works as expected.

Now, I would like to go further and be able to download in an automated way what can be obtained by clicking on one of the "download" buttons in the "Sequence" section of the above pages. More specifically, I'm interested in the "spliced + UTR" and "conceptual translation" fasta sequences.

If I use the DOM inspector in Firefox (Shift-Ctrl-C), I can locate things like the following:

    download

        >K06C4.12 spliced + UTR
        ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTCACCAAGCCAAAGGACGGAAAGAAGAG
        ACGTCATGCCCGTAAAGAATCCTACTCCGTGTACATCTACCGTGTGCTCAAGCAGGTTCATCCAGACACTGGAGTCTCTT
        CCAAGGCTATGTCCATCATGAACTCTTTCGTCAACGACGTCTTCGAGCGTATCGCTGCTGAGGCATCTCGTCTTGCTCAC
        TACAACAAGAGATCCACCATCTCTTCTCGCGAAATTCAAACTGCCGTCCGTTTGATTCTTCCAGGAGAGCTTGCCAAGCA
        CGCCGTCTCTGAGGGAACCAAGGCTGTCACCAAGTACACCTCTAGCAAGTAA

K06C4.12conceptual translation.fasta

        >conceptual translation
        MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH
        YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

    view conceptual translation (123 aa)

    >conceptual translation

    MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSSKAMSIMNSFVNDVFERIAAEASRLAH

    YNKRSTISSREIQTAVRLILPGELAKHAVSEGTKAVTKYTSSK*

I don't find these if I simply download the page with wget and search in the source code, so they must be dynamically generated.

How could I proceed to get the sequences?

The documentation about Wormbases's REST API is not very detailed, so I don't understand how to form an URL that would give me the desired sequences.

Further inspecting the source code of the downloaded page, I found the following:

    <span class="ui-icon ui-icon-close module-close" wname="sequences"

    Sequences

This suggested me to have a look at http://www.wormbase.org/rest/widget/transcript/K06C4.12/sequences and I obtained a page containing the elements found in the DOM inspector.

This can be parsed, for instance using the bs4 python module:

# To get the page
import requests
# To parse the page content
from bs4 import BeautifulSoup
# To get this into Biopython
from io import StringIO
from Bio import SeqIO

rest_template = "http://www.wormbase.org/rest/widget/transcript/{accession}/sequences"
rest = rest_template.format(accession="K06C4.12")
# This did not work with "html.parser"
soup = BeautifulSoup(requests.get(rest).text, "html5lib")

def generate_fasta(soup):
    name = None
    for div in soup.find_all("div"):
        if "ui-helper-hidden" in div["class"]:
            if div["id"] == "content":
                fasta = SeqIO.read(StringIO(div.text), "fasta")
                description_fields = fasta.description.split()
                if len(description_fields) == 4:
                    fasta.id = "_".join(description_fields[:2])
                    if name is None:
                        name = fasta.name
                elif len(description_fields) == 2:
                    # This assumes that the "conceptual_translation"
                    # is after the "unspliced + UTR" or "spliced + UTR"
                    # Otherwise, one would need to pass
                    # the accession as function parameter
                    fasta.id = f"{name}_translation"
                else:
                    fasta.id = "_".join(description_fields)
                yield (fasta.id, fasta)

records_dict = dict(generate_fasta(soup))
print(records_dict['K06C4.12_spliced'])
print(records_dict['K06C4.12_translation'])

Output:

ID: K06C4.12_spliced
Name: K06C4.12
Description: K06C4.12 spliced + UTR
Number of features: 0
Seq('ATGGCCCCACCAAAGCCGTCTGCTAAAGGAGCAAAGAAAGCCGCCAAGACCGTC...TAA', SingleLetterAlphabet())
ID: K06C4.12_translation
Name: conceptual
Description: conceptual translation
Number of features: 0
Seq('MAPPKPSAKGAKKAAKTVTKPKDGKKRRHARKESYSVYIYRVLKQVHPDTGVSS...SK*', SingleLetterAlphabet())


How to measure the total size of a fastq file in base pairs? Or Kbps/Gbps. It feels like it should be conceptually very simple, but I can't seem to figure out the right combination of keywords to find it via my search engine. Help would be appreciated!
I have BBMAP, SRAtoolkit and MEGAHIT already installed, and also use bash. I'd be very happy if this can be answered with software that I already have, but if not that's perfectly fine.

I've been using this:
cat file.fastq | paste - - - - | cut -f 2 | tr -d '\n' | wc -c
Explanation :

paste - - - - : print four consecutive lines in one row (tab delimited), to merge the info for each read

cut -f2 : print only the second column, to access the sequence after the paste

wc -c : count the characters

tr -d '\n': to remove from count the eventual newline characters

(a tip for your googling: try search for "counting number of bases in fastq file")


Demographic model for admixed African Americans I am fairly new to genomics and population genetics. And I am reading this article 'Demographic history and rare allele sharing among human populations':

https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3142009/pdf/pnas.201019276.pdf

In this article they propose a demographic model in Figure 4 on page 11986. However, I am confused how does African Americans fit in this demographic model that are known to have both African and European ancestry. And if I am to consider African American population than what is the split time between the two populations, time they merge into admixed population and the admixed proportion.
My apologies if it's a really basic question but, any insights would be appreciated.

Just to further my answer, the problem with molecular dating within a species (I don't do humans [excepting immunology], but it's the same theory), is it requires an assumption of clonality (bifurcating trees). Mendelian species recombine towards panmixia, which will make any attempt to date via de novo point mutation, appear far younger than they really are, or for panmixia, incalculable (it breaks the central assumption). Dating through admixture requires a lot of extra data and/or assumptions for example the generation time,  putative admixture per generation, and shifts in population dynamics. Maybe you can calculate some of those requirements but  alot will remain assumptions.

I have not read the paper you suggested, because I'm not involved in human genetics [outside immunology], but it is possible the authors may (and I stress may) have simplified an assumption that could be difficult to estimate. It is not clear what the objective of the paper was.

Quantitative analysis is extremely important in modern population genetics/ evolutionary theory because historically it was overlooked and resulted in lots of hypothese which were difficult to test.


How to get log2 fold change of RNA-Seq data for time series experiment? I know if there is one control and one treatment group it is pretty straight forward to interpret log 2 fold change. But, I have time course experiment. I have infected cells with viruses and I harvested cells at different time points: o hr, 6 hour, 12 hr and 24 hour. Here o hr is my non-infected control. What is the best way to show log 2 fold change for each gene in this type of time course experiment. I ultimate aim to find highly expressed genes at different time points.

If you just want the change with time then use time as a continuous variable. The log2FC will then be the change per hour (you can use the default Wald test to extract it).

For time-point specific comparisons you'll instead need to treat the time points as groups and use contrasts as appropriate.


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

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

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

Which is my error?
Thank you in advance.

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

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

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

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

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

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

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

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

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

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


how to generate an scRNA counts file I have barcodes.tsv, genes.tsv and matrix.mtx file from my RNA seq data and I want to generate a scRNA counts files using R. When I tried using the readMM function, it gave me the following error:

Error: cannot allocate vector of size 185.1 Gb.

This happens even when I run it as a job. How can I generate a combined csv file using R?

Thank you.

I use the DropletUtils package from Bioconductor. Getting counts begins like this:

library(DropletUtils)
filePath <- "path_to_data/raw_feature_bc_matrix/"
singleCellExperiment <- read10xCounts(filePath, col.names = TRUE)


isSpike function in SingleCellExperiment package is deprecated? I am trying to follow a tutorial from Sanger institute (from May 2019) on analysis of single cell RNA Seq data. They use isSpike function to filter out ERCC (control) and MT (mitochondrial RNA) reads, but I am getting an error message that isSpike is a deprecated function. I am trying to use altExp instead as recommended by some people. I can't seem to figure this out. Can anyone provide some code to use instead of isSpike? I am new to bioinformatics.

I tried to swap out the isSpike function with the altExp function, but I am getting stuck on how to translate the grepl function.

isSpike was provided by SingleCellExperiment. It is now deprecated. The notes in SingleCellExperiment indicate that one should use altExps

I am struggling with a similar problem. Can't compile a package (powsimR) because a dependency uses isSpike (not used in the code directly), so I have to hunt down where it is used by one of a very long list of dependencies.


Seurat define manually a cluster and find markers I am aware of this question Manually define clusters in Seurat and determine marker genes that is similar but I couldn't make tit work for my use case.
So I have a single cell experiments and the clustering id not great I have a small groups of 6 cells (I know it is extremely small, but nonetheless I would like to make the most of it) that are clearly isolate in UMAP and display marker that I am interresting in.
I managed to get those cells names in barcode (e.g.: AGCCAGCTCGCTTTAT).
I would like to compute differential expression of those cells against all others groups (using FindMarkers or FindAllMarkers).
I tried to manually define the cluster in:
Seurat_obj@meta.data$seurat_clusters
But using FindMarkers (it says this groups don't exist)
then I try to modify
s_19270BC@active.ident
but then it could not found metadata attached to my newly define group.
I am kind of stuck and any help is warmly welcome
Code
My R is quite rusty...
names_ # a vector with the names of the cell I want as a cluster

# my seurat object have 5 cluster 0..4. So I want to add a 6th clusters with id 5

#First attempt manually modify the vector active.ident and cluster
# they are both factor so we need to add a factor level
levels(seurat_obj@active.ident) = c(levels(seurat_obj@active.ident), 5)
levels(seurat_obj@meta.data$seurat_clusters) = c(levels(seurat_obj@meta.data$seurat_clusters), 5)
# then modify them :
seurat_obj@active.ident[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5
seurat_obj@meta.data$seurat_clusters[which(row.names(seurat_obj@meta.data) %in%  names_)] = 5

FindMarkers(s_19270BC, ident.1 = 5)

So This morning it works... I don't know what I have done yesterday...
EDIT:
I added the Code I used

It would have been better if you have shared your code as well as the error messages you got format these as such.

As far as I can understand, your problem is that you failed to make your cluster as default ids. For this you can use:

seurat_obj <- SetIdent(seurat_obj, id=seurat_obj@meta.data$column_of_interest).

The cluster identities are stored in the @meta.data slot of your object.


Why does SBOL have both Component and ComponentDefinition? What is the point of having both a Component and a ComponentDefinition in SBOL? Shouldn't you be able to link straight from one component to its sub-components?

The same question goes for Module and ModuleDefinition.

The key to understanding this is to think about a system with more than one copy of a genetic construct.

For example, let's say you have a construct that uses two copies of the pBAD promoter, each regulating a different coding sequence. Both of these promoters are pBAD promoters, with the same sequence and type information, so they should share the same ComponentDefinition.

But there are two copies. In the context of the larger construct, we need to be able to talk about each copy separately. For example, we might need to say that the first comes after a MoClo cut site, while the second comes after a terminator. For that, we need the Component object, to have a way to say "this particular copy of pBAD", as opposed to "pBAD in general", which is what the ComponentDefinition tells us about.


Get Gene Expression Matrix from GEOquery I am trying to get gene expression matrix for a list of genes
I have for my list: Gene ID, Gene Symbol
How can I get for each gene corresponding expression as array??
I suppose that I will have
GENE_ID GENE_NAME SAMPLE_1 SAMPLE_2 .....
I am new in that, I don't know more details
I am trying to study relation between diseases and genes(human) using gene expression as features.
Thanks in advance.


Is it OK to use Blast+ to query NCBI's 16s rRNA database with 16s DNA sequences? I have 16S DNA sequences from diversity sequencing, and want to query these against NCBI's 16S rRNA database.  If I ultimately want to get taxids for the species that my sequences have, does it make a difference that my sequences are 16S DNA and not RNA?  I am using Blast+.

I imagine this is how analysis of 16S sequences normally goes (get 16S DNA sequences, query them against a 16S DB), but does anyone know if Blast+ or the 16S rRNA DB itself allow for DNA queries?

In other words, if I were to blast a sequence of 16S DNA and its exact transcript RNA, would I get the same results?

Your confusion is understandable, since it isn't explained in the NCBI description of the database, but most of the sequences in the database are derived from DNA that encodes the 16S rRNA gene, not from the rRNA itself. In some cases, sequences are derived from cDNA that is synthesized using the actual rRNA as a template, but it's still DNA sequencing. Direct RNA sequencing is extremely uncommon. Depending on the cDNA sequencing protocol (first strand vs. second strand sequencing), the transcript sequence might be the reverse-complement of the 16S rDNA, but the only difference that you would see is that the numbering of the bases would be reversed. If you are sequencing rDNA from organelles, there may also be a few bases that are different, due to RNA editing.


Visualizing the difference among three groups I have 3 data frames for three groups of patients and in each of them I have the number of mutation types like insertion, deletion, SNP and total of mutations for each patients. In each group I have different number of patients like below

> head(dat1)
                patient  DEL  INS   SNP   total
1:    LP6008337-DNA_H06  927  773 40756   42456
2:    LP6008334-DNA_D02 1049  799 31009   32857
> dim(dat1)
[1] 21  6
>

 head(dat2)
                Patient   DEL  INS   SNP total
1:    LP6008031-DNA_E01 13552 3374 62105   79031
2:    LP6005500-DNA_G01   539  500 43451   44490
> dim(dat2)
[1] 33  6
>

> head(dat3)
                Patient   DEL   INS    SNP    total
1:    LP6005935-DNA_F03 39168 16739  58095   114002
2:    LP6008269-DNA_D08   849   910 103501  105271

> dim(dat3)
[1] 106   6

I want to show if the number of each mutational category and total is different between these three groups by chi-square test

Actual question is if total number of mutations, SNP, DEL and INS are statistically significant among groups. I used pairwise t test but I afraid this test is not a suitable test for the distribution of data moreover I don't know how to visualize the p-value

Imagining the comparison of total number of mutations between three groups, this picture is a good example

Or in this plot they compare several features but only two groups

Can you help ?

Thanks

Here a way to do it is to start first by creating a dataset containing your three different group, based on the few lines you display:

dat1$Condition = "dat1"
dat2$Condition = "dat2"
dat3$Condition = "dat3"
colnames(dat1)[1] = "Patient"
DF <- rbind(dat1,dat2,dat3)

             Patient   DEL   INS    SNP  total Condition
1: LP6008337-DNA_H06   927   773  40756  42456      dat1
2: LP6008334-DNA_D02  1049   799  31009  32857      dat1
3: LP6008031-DNA_E01 13552  3374  62105  79031      dat2
4: LP6005500-DNA_G01   539   500  43451  44490      dat2
5: LP6005935-DNA_F03 39168 16739  58095 114002      dat3
6: LP6008269-DNA_D08   849   910 103501 105271      dat3

Then, you can reshape DF in order to put all values DEL, INS, SNP and Total into a single column and create a categorical column. Doing this will be helpful for later producing the facets of your plot. To do that, you can use pivot_longer from tidyr package:

library(tidyr)
DF <- DF %>% pivot_longer(.,cols = c(DEL,INS,SNP, total), names_to = "var", values_to = "val")

# A tibble: 24 x 4
   Patient           Condition var     val
   <chr>             <chr>     <chr> <int>
 1 LP6008337-DNA_H06 dat1      DEL     927
 2 LP6008337-DNA_H06 dat1      INS     773
 3 LP6008337-DNA_H06 dat1      SNP   40756
 4 LP6008337-DNA_H06 dat1      total 42456
 5 LP6008334-DNA_D02 dat1      DEL    1049
 6 LP6008334-DNA_D02 dat1      INS     799
 7 LP6008334-DNA_D02 dat1      SNP   31009
 8 LP6008334-DNA_D02 dat1      total 32857
 9 LP6008031-DNA_E01 dat2      DEL   13552
10 LP6008031-DNA_E01 dat2      INS    3374
# … with 14 more rows

Now, you can plot your data and use the function stat_compare_means from ggpubr package to get the significance display on each tab. Based on the size of your group, you may want to adjust for the type of statistical test to be used (take a look at ?stat_compare_means to find the statistical test adequate to your dataset). You can pass various conditions to test as a list by using comparisons argument.

For facetting your graph, you can use facet_grid (orfacet_wrap, try the one you preferred). Altogether, you can get the plot by writing:

library(ggplot2)
library(ggpubr)
ggplot(DF, aes(x = Condition, y = val, fill = Condition))+
  geom_boxplot()+
  stat_compare_means(comparisons = list(c("dat1","dat2"), c("dat1","dat3"), c("dat2","dat3")))+
  stat_compare_means()+
  facet_grid(var~., scales = "free")

Hope it helps you to figure it out how to represent your data.


Aligning many long sequences I'm faced with having to align many (some 100s) bacterial genomes, where the genome length is in the millions. Obviously, this is beyond normal alignment techniques and it's unclear to me what the best practice is for such circumstances:

conventional alignment on a very powerful computer with lots of memory
break up the genome into smaller fragments and align them individually
some exotic different procedure

What possible avenues of attack are there?

(I've attempted to use Mafft and Clustal with little success)

Whole genome aliment can be done using Progressive Mauve, LAST or Mummer. For bacteria I used Mauve since it has also very nice visualisation engine. A very new tool is Minimap2, a super fast mapper that supposed beside read mapping be able to handle reference vs reference. However, I do not know how performance of it compares to the tools mentioned above.

If you are interested in a rough idea of the shared genome regions, you can use bevel. Bevel is not really an aligner, it is more like a dot-plot, but it is super fast (even for mammalian sized genomes).


Python string editing in Snakemake I'm running a custom perl script using Snakemake. I use this rule:

rule complexity_20mer_counter:
    input:
       os.path.join(fastq_trimmed_dir, '{sample}_{read}' + fastq_extension.split('.')[0] + '_val_' + '{read}'.strip('R') + '.fq.gz')
    output:
        os.path.join(custom_qc_dir, '{sample}_{read}_20mer_counter.txt')
    shell:
        '''
        #!/bin/bash
        perl scripts/20mer_counter.pl {input} > {output}
        '''

It fails (I'll get to why later) with this error:

Building DAG of jobs...
MissingInputException in line 1 of /rst1/2017-0205_illuminaseq/scratch/swo-406/2019-01-24/rules/20mer_counter.smk:
Missing input files for rule complexity_20mer_counter:
fastq_trimmed/xs0067_P2018SEQE27R01_S1_R1_001_val_R1.fq.gz

... this is because I use the output of trim_galore, this rule looks like this:

rule trim_galore:
    input:
       [os.path.join(fastq_raw_dir, '{sample}_'+ read + fastq_extension)
       for read in reads]
    output:
        timmed_fastq_files = [os.path.join(fastq_trimmed_dir, '{sample}_' + read + fastq_extension.split('.')[0] + '_val_' + read.strip('R') + '.fq.gz')
        for read in reads],
        trimming_reports = [os.path.join(trim_galore_qc_dir, '{sample}_'+ read + fastq_extension + '_trimming_report.txt')
        for read in reads]
    conda:
        "../envs/fastqc.yaml"
    shell:
        '''
        #!/bin/bash
        trim_galore --fastqc --gzip -o {fastq_trimmed_dir} --paired --retain_unpaired {input}
        # Move all qc reports to the fastq_trimmed_dir
        mv {fastq_trimmed_dir}/*.zip \
           {fastq_trimmed_dir}/*.html \
           {fastq_trimmed_dir}/*.txt \
           {trim_galore_qc_dir}
        '''

As you can see, trim_galore not only uses my {read} parameter in the output filename but again adds it at the end, only now without the "R".

I tried solving it in my first rule by applying a python .strip('R') on my {read} wildcard, but this does not seem to work. How can I correctly get '1' and '2' from 'R1' and 'R2' (which are my {read} wildcards)?

(Extra bonus points for a nice output.trim_galore.something as input for my complexity_20mer_counter rule, which does not work because the output is a list (it needs 2 fastq files at the same time), but my complexity_20mer_counter rule can be run in parallel again.)

It seems that it is not possible to apply strip() directly on the wildcard placeholders. Luckily it is possible to executed functions in the input directive.

So try changing your input of the complexity_20mer_counter rule to this:

lambda wc: os.path.join(fastq_trimmed_dir, wc.sample + '_' + wc.read + fastq_extension.split('.')[0] + '_val_' + wc.read.strip('R') + '.fq.gz')

The wildcard object is passed to our anonymous function. So have access to the wildcards value via wc.sample and wc.read, which are now handled as normal strings and we can apply strip on them.


Can I index a compressed FASTA file using STAR? I am using STAR to align RNA-seq reads to a reference genome. Before the alignment, I need to generate an index of the reference genome. I use the following code to generate the index successfully:

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa --sjdbGTFfile ref.gtf --sjdbOverhang 100

This works fine. However, I would like to keep my reference genome compressed to save disk space. So I am trying the following command:

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star --genomeFastaFiles ref.fa.gz --readFilesCommand "gunzip -c" --sjdbGTFfile ref.gtf --sjdbOverhang 100

but I get the following error:

EXITING because of INPUT ERROR: the file format of the genomeFastaFile: ref.fa.gz is not fasta: the first character is '' (31), not '>'.
 Solution: check formatting of the fasta file. Make sure the file is uncompressed (unzipped).

I am using the readFilesCommand successfully with compressed RNA-seq fastq files. Does anybody know if there is a similar option to use compressed references? Is there a workaround using Unix commands (maybe piping?) or do I need to decompress the reference, index it and then compress it again?

I can think of two possible workarounds. The simplest would be to try using zcat instead of gunzip -c, in case that works. Assuming it doesn't, you can use process substitution, if your shell supports it. Something like this:

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles <(zcat ref.fa.gz) --sjdbGTFfile ref.gtf --sjdbOverhang 100

The <() creates a file descriptor and passes that to the command. It is very useful for passing streams to commands that expect files.

If your shell doesn't support this construct, assuming you are on a system that supports named pipes (Linux does), you can do the same thing a bit more manually:

Create a named pipe (FIFO)

mkfifo foo

Set it to read the uncompressed genome

zcat ref.fa.gz > foo &

The & is important since the command above will hang until something starts reading from foo. Alternatively, you can just run this in one terminal and the next command in another.
Read the input from the pipe

STAR --runThreadN 8 --runMode genomeGenerate --genomeDir output/index/star \
    --genomeFastaFiles foo --sjdbGTFfile ref.gtf --sjdbOverhang 100


Why the t-test for a specific gene shows different value compared to differential analysis? I have RNA-Seq data for LUNG cancer. 370 tumor and 50 Normal. For differential analysis initially I did some filtering and kept approx. 19k genes for further analysis. I used edgeR.

With a FC > or < 2 and FDR < 0.05
tr <- glmTreat(fit, contrast=contrast.matrix, lfc = log2(2))
I got approx 3000 differentially expressed genes.

Among the 3000 differentially expressed genes I'm interested in a gene OIT3.

GeneSymbol  logFC    unshrunk.logFC   logCPM          PValue      FDR
  OIT3  2.156487448  2.171599487    -0.437234244    0.000701037 0.005304312

I see the above gene OIT3 is Upregulated. Why I see logCPM with negative value?

And I plotted the logCPM expression of the above gene and did a t-test which gave pvalue 0.1 which means not significant.

1) I see the gene is differentially expressed with Differential analysis using edgeR. Why I see logCPM with negative value?

2) With t-test why it is not significant?

There is no reason your t-test should reproduce edgeR. In fact, edgeR exists because t-test is inappropriate.

edgeR does the tests by pooling information from all genes, because with the low number of replicates your t-test doesn't give sufficient statistical power.

What you need to do is: check visually your gene and make a decision yourself based on graphical results and edgeR results if your gene shows biologically significant.


Doubt about using TPM for statistics We designed an experiment to explore the potential role of carbon dioxide on algae physiology using RNA-Seq. We analyse the differential gene expression using DESeq2 but now we are interested into analyse a few genes in detail.
My PI say to make the statistics over the ration response (expression_treatment/expression_control - 1) and my question is the following: Can I use TPM to calculate this ratio and make following statistics? I read in the documentation of DESeq2 that TPM we can't use this normalisation method to differential expression and my my other question is Why shouldn't?
My advisor ask to re-do the statistical analysis using the ratio response (rr) so i need the rr for each sample. I think this doesnt make sense cause the result will be or should be the same that using DESeq2.
Update: I'm using the fold-change calculated by DESeq2 and I'm not sure if I have transformed the data correctly. I took the two columns log2FoldChange and lfcSE (both are log-transformed variables) so I went back to the non-transformed variable as follow: FC=2^log2FoldChange and SE=lfcSE*2^log2FoldChange. I saw this approximation in this post. Is that correct? I have obtained a very broad range of error/

This is a common question, and the answer is while you can calculate this, it's not statistically robust.
You're likely to arrive at false conclusions.
Why is this?
Because while TPM is something you can calculate from a single sample, your question involves multiple replicates and those replicates vary between each other.
Counts in units of TPM are missing information about the variation between samples.
This extra information affects fold change estimates, and is why differential analysis of sequencing data is complicated.
See this answer to a related question, and this often-cited blog post for more information.


Bash scripting FastQC for multiple fastq files in multiple directories I am completely new to bioinformatics so I'm looking to learn how to do this.

I have multiple directories with fastq files: E.g; 10 Directories with each time series, each with Treatment and control directories, each with rep1 rep2 rep3.

For example: T9/Infected/Rep1/*.fastq.gz.

I'm looking to create a loop to run fastQC on each fastq file instead of having to submit a separate job for each directory.

Then to either output the fastQC data to a single directory or if possible a directory corresponding to each rep - e.g. rep1 results go into a folder called rep1 and so on.

multiqc kind of glazes over some important information, like the exact adapters and duplicated sequences in a library. If you plan to spend big $$ for sequencing a library it is better to look at both the multiqc report and the actual fastqc html report to get a better idea of any error modes.

Going off of @Kubator's answer, I noticed that there was no command to run fastqc.

Here's a simple one-liner to run fastqc in parallel on all of your fastq files. The -j 25 uses 25 threads. Change 25 to however many threads you want/have for max speed.

#Run fastqc on everything in parallel.
> find ../reads/ -name '*.fastq.gz' | awk '{printf("fastqc \"%s\"\n", $0)}' | parallel -j 25 --verbose
# copies all the fastqc files to directory ./
> find ../reads/ -name '*fastqc.*' | xargs -I '{}' mv '{}' ./

These files might be output by multiqc, anyway!

edit - now the command is safe for whitespace in file names. Please, no newlines in filenames though!


Best practices for deciding if two structural variants are actually the same variant? I know that I can use bcftools isec to compare two VCF files containing single nucleotide variants (SNVs); it will generate four possible output files: one with all of the unique called SNVs in file 1, one with all of the unique called SNVs in file 2, one with all of the SNVs from file 1 that were also found in file 2, and finally one with all of the SNVs from file 2 that were also found in file 1.

Is there a similar type of tool that will allow one to compare structural variants (SVs) between two VCF files?  I am about to tweak one of my Perl scripts (that currently uses the VCF.pm Perl module to parse VCF files) to perform this calculation (because thus far I have not found any tools designed for SVs) but I thought I would ask in this forum in case anyone else has already invented this wheel.

The tool that you want to use for this task is bedtools window.


removing nextera transposase adapters, cutadapt I am learning about NGS analysis and im currently learning about QCing and removing adaptors.

I am working on SRR1972920_1.fastq file.

When running fastqc tool on that file, adapter contamination is present in the form of Nextera Transposase adapters.

Searching for that adapter sequence via google the sequence for at the 3' end to be removed is TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG (obtained via illumina adapter sequence-Nature).

Running a very basic cutadapt command line to remove the adapter

cutadapt -a TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG -o <any output name> <my input file>

After running the fastqc tool again on the new output file from cutadapt, it seems nothing changed and the contamination still present.

However, when running the same command above but changing the sequence to CTGTCTCTTATA, which is the sequence of transposase adaptors according to fastqc tool adapters_list.txt file all is good after running fastqc tool on the file

cutadapt -a CTGTCTCTTATA  -o <any output name> <my input file>

Now, why this is happening? which Nextera transposase sequence is the correct? or am I doing something wrong?

The adapter sequence you have googled is the sequence on the adapter primer. This works when you want to remove primer-dimers.

With transposase adapters or ATAC seq, you have very short fragments and you sequence into the adapter. So what you is sequenced is actually the reverse complement. See below image from trimmomatic manual

So if you look at what worked, and the sequence you googled, they are reverse complement of one another.

You can also see this in the NexteraPE-PE.fa, Trans1 and Trans1_rc @Haci recommended:

more NexteraPE-PE.fa
>PrefixNX/1
AGATGTGTATAAGAGACAG
>PrefixNX/2
AGATGTGTATAAGAGACAG
>Trans1
TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG
>Trans1_rc
CTGTCTCTTATACACATCTGACGCTGCCGACGA
>Trans2
GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG
>Trans2_rc
CTGTCTCTTATACACATCTCCGAGCCCACGAGAC

To answer your question, if most of your fragments are short, and there's not a lot of primer dimers, then yes you are doing the right thing.


Difference between paired-end, mate-pair and long read I writing here because I have some questions for you.
I wondered what the essential differences were between paired-end, mate-pair and long read?
for me mate-pair and classic paired-end are both paired-end reads, with the difference that :

For classical paired-end: the insert size on classic paired-end is
smaller (about 500bp)
while the insert size of mate-pair is much longer (several Kb)
which allows to join the contiguous between them especially is it?

In revenge for the long-reads, I imagine that they are simply reads that are synthesized with a large read size but that do not allow like the maite-pair to make junctions between contigs?

They are all very different in separate regards, but they all refer to different wet-lab and sequencing protocols/technologies.
First, PE (paired end) reads are typically short (50-300) reads, most often Illumina HiSeq, MiSeq or NovaSeq protocols. Both pairs originate from a single fragment which is sequenced from either end:
------> PE1
-------------------------------- Fragment
                         <------ PE2
       ^^^^^^^^^^^^^^^^^^        Insert

In contrast mate pairs arise from a fragment that is circularized before sequencing. Now imagine the fragment is circularized:
                     -----> MP1
....--------------------------------.......
          MP2 <-----
                    ^Circularized (joined) here

Long reads typically refer to PacBio or Oxford nanopore data, with very long reads lengths (up to several hundred kb).


Get gene sequence based on the annotation I've got the reference genome with Python like so:
sequences_by_chr = {}

with gzip.open("data/Zea_mays.B73_RefGen_v4.dna.toplevel.fa.gz", "rt") as f:
    for seq in SeqIO.parse(f, "fasta"):
        sequences_by_chr[seq.id] = seq.seq

I've also parsed a GFF3 annotation and would like to get the gene sequence based on said annotation.
Suppose the following:
chromosome = 1
strand = "-"
start = 45286331
end = 45290076

Since it's on the "minus" strand, should I do the following:
sequences_by_chr[chromosome].complement()[start:end]
or should I use the reverse_complement()?
I'm completely confused by this, would appreciate any hint!

You can parse the GFF3 to generate a bed-like file with gene ranges and use that with bedtools getfasta to obtain a multi-FASTA file with gene sequences. An example case is shown below:
## get fasta and gff3 files
wget ftp://ftp.ensembl.org//pub/current_gff3/saccharomyces_cerevisiae/Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz
wget ftp://ftp.ensembl.org:/pub/current_fasta/saccharomyces_cerevisiae/dna/Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## unzip fasta
gunzip Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa.gz

## generate genes.bed
zgrep -v '^#' Saccharomyces_cerevisiae.R64-1-1.100.gff3.gz | awk 'BEGIN{FS="\t";OFS="\t"}($3=="gene"){print $1,$4-1,$5,"name",1000,$7}' > genes.bed

## create gene fasta
bedtools getfasta -fi Saccharomyces_cerevisiae.R64-1-1.dna.toplevel.fa -bed genes.bed -s -fo genes.fa
```


Diploid Consensus Sequences My goal is to generate a diploid consensus genome. By "diploid consensus", I mean that I want to merge all supporting reads for a genomic region into a single haplome such that the resulting sam file contains TWO sets of reads, one from each haplome.
I have Illumina short-read files from the NCBI's SRA. These short-reads were obtained from a diploid species.
I am worried that, by simply using Picard's MarkDuplicates with the "REMOVE_DUPLICATES" option set to "TRUE", I will retain only the homolog with the highest quality and output a haploid version of my genome.
Here is my current workflow:

## Map Illumina reads to the reference genome. Include read IDs
bwa mem -M -R '@RG\tID:BioSample.LibraryName\tSM:BioSample\tLB:LibraryName\tPU:SingleFlowcell.SingleLane\tPL:ILLUMINA\t' reference.fasta SRR[number]_1.fastq SRR[number]_2.fastq > SRR[number].sam

## Sort Illumina reads by coordinate
java -jar $EBROOTPICARD/picard.jar SortSam I=SRR[number].sam O=sorted_SRR[number].sam CREATE_INDEX=true SORT_ORDER=coordinate TMP_DIR=$TMPDIR

## Remove duplicate reads from sam file
## *Note: How do I remove sampling error reads?
java -jar $EBROOTPICARD/picard.jar MarkDuplicates REMOVE_DUPLICATES=true REMOVE_SEQUENCING_DUPLICATES=true I= 'sorted_SRR[number].sam' O= nodup_SRR[number].sam METRICS_FILE= metrics_SRR[number].txt TMP_DIR=$TMPDIR

## Convert duplicate-free sam file to bam file
samtools view -S -b nodup_SRR[number].sam > nodup_SRR[number].bam

## Re-align reads around indels to reduce slippage.
gatk LeftAlignIndels -R reference.fasta -I nodup_SRR[number].bam -O aligned_SRR[number].bam

Does anyone know of a program which will whittle my sam files down to ONLY two homologs per region?

How to separate reads, one for each reference haplotype.

Map all the reads against the reference genome
Use GATK to call SNPs
Use GATK tool [FastaAlternateReferenceMaker][1] to build the alternative reference
You now have the two "haploid references" (quotes needed because you have to take into account that they will not necessarily reflect the exact haplotypic structure, it could be a mix of the two haplotypes)
Align the reads against the two references, and attribute each read to the reference to which they map better. My colleagues developed for a different purpose (detection of Allelic Imbalance) a pipeline that includes the step for determining to which of two references a read belongs. The tool is called [BayesASE][2]

In general you may study all the work that has been done in the field of allelic imbalance and allele specific expression, where this problem has been encountered and solved.
[1]: https://gatk.broadinstitute.org/hc/en-us/articles/360037594571-FastaAlternateReferenceMaker
[2]: https://github.com/McIntyre-Lab/BayesASE