How to merge .fastq.qz files into a single .fastq.gz with their same id without losing any content in parallel

All you need is cat. You won't find any better tool for a simple job like this. Just run:

cat NA24694_GCCAAT_L001_R1*fastq.gz > EA00694_GCCAAT_L001_R1.fastq.gz
cat NA24694_GCCAAT_L001_R2*fastq.gz > EA00694_GCCAAT_L001_R2.fastq.gz
cat NA24694_GCCAAT_L001_R3*fastq.gz > EA00694_GCCAAT_L002_R1.fastq.gz
cat NA24694_GCCAAT_L001_R4*fastq.gz > EA00694_GCCAAT_L002_R2.fastq.gz

That should work just fine. You say you tried cat and it "didn't work" but since you don't tell us how it failed, I can't really help. The only issues I can think of is that either there are too many files (as in several hundred thousand, whatever the value returned by getconf ARG_MAX on your system which is 2097152 on mine) or, more likely, you are running out of disk space.

If it's a disk space issue, you might be able to get around it by adding each file and then deleting it:

for i in 1 2 3 4; do
    for file in NA24694_GCCAAT_L001_R${i}_*fastq.gz; do
        cat "$file" >> EA00694_GCCAAT_L001_R${i}.fastq.gz && rm "$file"
    done
done

IMPORTANT: the command above is destructive. It will delete each file after it has been added to the new one. If everything works fine that's not a problem, but if I've made a mistake or your file names are slightly different you might lose your data. So I strongly urge you to simply run the cat commands in the first section on a machine with enough disk space to deal with them.


How does htslib/samtools access optional BAM fields?

How does htslib access each of the optional fields, given this structure?

htslib uses a macro, bam_get_aux, for this purpose. You can check the source code of this macro — it essentially just advances the data pointer of the current record to the start of the optional data, which is after all other data. But you shouldn’t care about this, and just use the relevant exported functions for reading optional (“auxiliary”) data:

bam_aux_get returns a pointer to an optional record, given by its two-character tag name.
bam_aux2‹C› converts the pointer returned by the previous function to the actual type of the optional tag; here, ‹C› is a type identifier: i for an integer, Z for a zero-terminated character string, etc. Refer to the documentation and the htslib/sam.h file for a complete list.

  Does it keep the order of the optional fields, or is this information unavailable?

As long as you don’t modify optional fields, the order is kept, but unspecified. Do not rely on the order!

  When a user runs htslib or even samtools view, how does the BAM row "end"

BAM records store a header that contains the record size.


What are the ways to keep track of branches in the analysis?

Save the different scripts with git (seems overkill)

Whoa. I did an actual double take when reading this:1 it’s the opposite of overkill.

Version controlling your scripts (using Git or something similar) is the absolute minimum, and should become completely automatic. For every new project I begin, one of the very first steps is to issue the git init command, and to set up a remote repository (on Github).

To keep track of different analyses, I use a combination of the following approaches:

Write reusable functions/scripts and parametrise. The parameters are kept either inside the script itself (that then calls the relevant function repeatedly), or in a Makefile (I recommend Snakemake).
Document the alternative analysis approaches; once again, this could be a Makefile with different rules for alternative analyses, or a set of notebooks (via R Markdown).
Have different Git branches for mutually exclusive approaches. At the end of the analysis one of these branches gets merged into master, and published. If I want to publish several analysis approaches, I merge all these branches into master, and use approaches (1) or (2) enable them simultaneously.

In fact, I recommend creating a Makefile for every analysis; I have found that this is the most practical, self-documenting way to run an analysis. It most closely resembles a wet-lab lab notebook. The advantage over a single R Markdown document is that rerunning just parts of the analysis can be completely automated, and dependencies in the workflow are apparent from the dependencies of the Makefile rules. This is much harder in R Markdown.

Some time ago I create an example analysis workflow to show how this can be structured. Nowadays I would use Snakemake instead of GNU make.

Regarding your other point:

  Make notes in the script itself

“Notes” are a dangerous beasts: documentation is important, but experience shows that it’s sometimes very hard to keep documentation synchronised with the code. There is no mechanism that ensures that documentation and code actually agree. Differences between presumed analysis and actually executed analysis can become very problematic.

People therefore prefer using self-documenting code as much as possible; that is: writing code so that its meaning becomes immediately clear, without comments, even to somebody who hasn’t worked on the code before. Doing this well is hard and takes practice, but improves overall code quality. Once again, using a Makefile helps here because the dependencies between the rules are self-documenting the kind of analysis that was performed.

Jeff Atwood has written two seminal essays on this subject:

Coding without comments, and
Code tells you how, comments tell you why.

They are two of the best pieces of advice on programming that I could give.

1 To emphasise: take a look at the edit history of this answer.


this is the swift standard primer file. You are a bit off in your interpretation of the columns (seems like you missed one as there are 12 columns, not 11 which threw everything off):

1) Chromosome

2) Target Start

3) Target Stop

4) Target Name (change which affected everything downstream in counting)

5) 5' Primer Start

6) 5' Primer Stop

7) 5' Primer Name

8) 3' Primer Start

9) 3' Primer Stop

10) 3' Primer Name

11) 5' Primer Sequence

12) 3' Primer Sequence

To convert to a BEDPE describing the primer locations, you only will need columns 1, 5, 6, 8, and 9 (and optionally 4). The last letter of the Region name for the denotes the strand orientation of the region. For example Region 1 has _F denoting forward or (+) strand for BEDPE usage, while region 2 has _R in name denoting reverse or (-) strand for BEDPE. Putting all this together, you can convert to BEDPE with:

$ awk -v OFS="\t" '{
            o=".";
            if($4~/_F\//){
                o="+"
            }
            else if($4~/_R\//){
                o="-"
            }
            print $1,$5,$6,$1,$8,$9,$4,".",o}' foo.masterfile
1   43814932    43814951    1   43815064    43815081    MPL_S505N-A519T_F/2 .   +
1   115256457   115256480   1   115256576   115256602   NRAS_G60-S65_new_R/1    .   -


Seurat has a vignette specifically for combining multiple 10x libraries.


You will probably be interested in the following UCSC wiki page, which explains how to go from most of the UCSC tables to GTF/GFF:
http://genomewiki.ucsc.edu/index.php/Genes_in_gtf_or_gff_format

The basic gist is that UCSC doesn't store any data internally as GTF or GFF, and so you will need to use our genePredToGtf utility in order to convert from our internally stored genePred format into GTF/GFF3.

This is the command I usually give when we get this question on our mailing list:

$ mysql --user=genome --host=genome-mysql.cse.ucsc.edu -A -N -e "select * from refGene" hg19 | cut -f2- | genePredToGtf -source=hg19.refGene.ucsc file stdin stdout

Change stdout to the output filename you want in the last command to get an hg19 refGene GTF file:

chr1    hg19.refGene.ucsc   transcript  11869   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357";  gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    11869   12227   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "1"; exon_id "NR_148357.1"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    12613   12721   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "2"; exon_id "NR_148357.2"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   exon    13221   14362   .   +   .   gene_id "LOC102725121"; transcript_id "NR_148357"; exon_number "3"; exon_id "NR_148357.3"; gene_name "LOC102725121";
chr1    hg19.refGene.ucsc   transcript  11874   14409   .   +   .   gene_id "DDX11L1"; transcript_id "NR_046018";  gene_name "DDX11L1";
...

If you have further questions about UCSC data or tools feel free to send your question to one of the below mailing lists:

General questions: genome@soe.ucsc.edu
Questions involving private data: genome-www@soe.ucsc.edu
Questions involving mirror sites: genome-mirror@ose.ucsc.edu


In case an answer would still be relevant after a year and a half, below is a solution with the Biostrings package.

library(Biostrings)

# loading data.table for the fread() function.
# fread() is not only fast but is quite good at guessing how
# the file is formatted
library(data.table)

# readDNAStringSet() function creates a "DNAStringSet" object
my_fasta <- readDNAStringSet("fasta.fa")

> my_fasta
  A DNAStringSet instance of length 4
    width seq                                                                names
[1]    61 TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG      rnd-5_family-5445...
[2]   485 AATCTGAGATTAGATTTATTTATCTGTTATCC...GGCCAATCATTCAAAAGTTATGAAAGTGATC rnd-5_family-2614...
[3]   215 TTTCTAGTTTGTATTTGAGGAAAGAAAATCCT...TAGAAATTGTAATTTACAAAATGAATTACTT rnd-5_family-2614...
[4]    87 ATGCTTAATTTGCTAAGCATTTAATAACCTTT...GTAGACATTAAAATAAACCTTATACTGCTAT rnd-4_family-798_...

# the fasta header is stored at the @NAMES slot of the @ranges slot of the
# DNAStringSet object created above
> my_fasta@ranges@NAMES
[1] "rnd-5_family-5445_Unspecified(-)" "rnd-5_family-2614_Unspecified(-)"
[3] "rnd-5_family-2614_Unspecified(-)" "rnd-4_family-798_Unspecified(-)"

# removing the trailing "_Unpecified..."
my_fasta@ranges@NAMES <- sub("_Unspecified.*",
                              "",
                              my_fasta@ranges@NAMES)

my_coord <- fread("repeat_masker.txt")
> my_coord
           V1   V2   V3                            V4 V5 V6
1: KB824701.1  417  478 rnd-5_family-5445_Unspecified  .  -
2: KB824701.1  587 1072 rnd-5_family-2614_Unspecified  .  -
3: KB824701.1  914 1129 rnd-5_family-2614_Unspecified  .  -
4: KB824701.1 1138 1225  rnd-4_family-798_Unspecified  .  -
#removing the ">" sign

my_coord$V1 <- sub(">",
                   "",
                   my_coord$V1)

# appending first 3 columns of the coordinate file to "NAMES" as requested by the OP
my_fasta@ranges@NAMES <- paste(my_fasta@ranges@NAMES,
                               my_coord$V1,
                               my_coord$V2,
                               my_coord$V3,
                               sep = "_")

# writing the modified DNAStringSet object into a fasta file
writeXStringSet(my_fasta, "edited_fasta.fa")

>rnd-5_family-5445_KB824701.1_417_478
TTTCACCGTAAATTAACTTTGAGAGGAGCTAATTCCTAAAAGAATTATACCGGCCTATTTG
>rnd-5_family-2614_KB824701.1_587_1072
AATCTGAGATTAGATTTATTTATCTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTA
CCACTTTTGATTCATTTTTTAAATCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTTGA
AGGCAACACGAATCCAGAGTGATGCTTTACATAAATCTGCTTCTACCGATGCCAAAAATTGACGATATTCTATTATTTAA
TCTAAATGTTAGTCTTTACATACCCTCCCCTAATTGTTAGAATTTTATGAAATTTGATTTCAGGGGTCAGTTTAGCATGC
TAAATCTAATTCAATAGATTGATATTTTTCTTCAGGTAAAGAAAATTTTTGCGTCAAAGTAATCATATTCCTCCACGATT
GCATATAACTATGGTATATAATTTAAAAGATTACACTTTACGTAATGAAAAATCGGCCAATCATTCAAAAGTTATGAAAG
TGATC
>rnd-5_family-2614_KB824701.1_914_1129
TTTCTAGTTTGTATTTGAGGAAAGAAAATCCTACTATCACTTTTTAAAAAATTCACCAATCTGAGATTAGATTTATTTAT
CTGTTATCCGCTCAGGTTGAAAAGTTTGCGCAATGTATAACTATAAAAATCTTTCTACCACTTTTGATTCATTTTTTAAA
TCTGGGGTCACATTTTATTCACGGTAGAAATTGTAATTTACAAAATGAATTACTT
>rnd-4_family-798_KB824701.1_1138_1225
ATGCTTAATTTGCTAAGCATTTAATAACCTTTTTGGGTTTTGTGATAAAGGATGCTGTAGACATTAAAATAAACCTTATA
CTGCTAT


Why does Coronavirus genome contain poly A region at the 3' end in GenBank?

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


R: 'Matrix' can not be unloaded, but 'writeMM' method not found

Ultimately I just restarted Rstudio as Konrad Rudolph suggested, imported Matrix library, and then writeMM method worked just fine. However, rownames were converted into just numbers, but I needed to keep the gene names instead as it was in dgcMatrix, so after some search I found a different way of doing that from here:

http://www.sthda.com/english/wiki/saving-data-into-r-data-format-rds-and-rdata

So, I just did save() (and then, to load back - load()) and the rownames are maintained like I wanted it to be. So, storing it now as .RData file.


How to download multiple proteomes at once?

You could try download-refseq-genomes. It will fetch all genomes from the NCBI FTP server that are in a specific subtree of the phylogeny.

For example, downloading the amino acid sequences from all proteins in each genome assembly for all species belonging to Chloroflexi (NCBI taxon ID = 200795):

download-refseq-genomes.pl -t faa 200795

This will download all assemblies with status=Complete. Use option -a to download all assemblies with any status.


You have an obvious error in your shell script that should have given you an error message:
sra_code = "DRR163""$i"

Should be:
sra_code="DRR163""$i"

You cannot have spaces around the = sign in a variable assignment. The line sra_code = "DRR163""$i" means "run the command sra_code with the arguments = "DRR163""$i". Since there is no such command, you see an error.
Here's a working version of your command:
for i in {256..345}; do
    sra_code="DRR163$i"
    fasterq-dump "$sra_code" -p -v --outdir ./output/ -e 8
done


Is there a Python/R package with the ability to convert an alignment and reference into a CIGAR?

Your request is easy to implement. I wouldn't use any libraries in this case. As you haven't showed your implementation, I will provide one:

def gen_cigar(ref, qry):
    if len(ref) != len(qry):
        raise Exception('unequal length')
    cigar = []
    for i in range(len(ref)):
        r, q = ref[i], qry[i];
        if r is '-' and q is '-':
            raise Exception('both gaps')
        op = '=' if r is q else 'I' if r is '-' else 'D' if q is '-' else 'X';
        if len(cigar) > 0 and cigar[-1][1] is op: # add to the last operation
            cigar[-1][0] += 1
        else: cigar.append([1, op]);              # a new operation
    return "".join(map(lambda x: str(x[0]) + x[1], cigar)); # turn to string

print(gen_cigar('ACGTAT-CT', 'ACGT-TGGA'));

The output is:

4=1D1=1I2X

PS: It is interesting that you often approach problems from a different angle. The solutions that you think "must exist" actually rarely exist because those are not the ways we (me at least) work with data. This is why I prefer an asker to explain the bigger picture in addition to isolated intermediate problems – there could be simpler and more straightforward solutions without the need of solving the contrived problems.


Using this data frame:

gaps
  Chrom Start      End    Name
1     I     0 15072434 contigs
2    II     0 15279421 contigs
3   III     0 13783801 contigs
4    IV     0 17493829 contigs
5     V     0 20924180 contigs
6     X     0 17718942 contigs
7 MtDNA     0    13794 contigs

and this command:

chromPlot(bands = gaps)

I could obtain this plot:

and then, if you want to add some segments for a subset of bases, you can do:

gaps2 <- gaps[1:2, ]
gaps2$End[2] <- 7000000

chromPlot(bands = gaps, segment = gaps2)

and obtain the following plot:

A more satisfactory answer could put segments of different colour on the chromosome, so bands but this is the best I've got so far.


Galaxy has API and API-consuming libraries (such as BioBlend) that will allow you to interactively script against it without opening the graphical interface at all.

However you can also take almost any tool out of Galaxy and use it independently since everything is open source. The converter you mentioned is available as a Python script here and the tool 'wrapper' which you can use to understand how to invoke the python script is next to it.


How to find all variable-length seqs with an exact 5' and 3' match in a FASTA file

python 3 solution

This python program will get you what I think you want. It is not elegant given that it (A) stores the whole sequence in memory and (B) spends most of its time slicing strings rather than working with a buffer. This program is not recommended for FASTA files for which the largest sequence will not fit well in RAM.

All that you need to do to modify it for other genomes or other binding sites is to change the three lines that begin with forseq, revseq, and filename.

from Bio import SeqIO
from Bio.Seq import Seq

forseq = Seq("tttaca".upper())
revseq = Seq("gtcttt".upper())
filename = "NC_003030.fasta"

forcom = forseq.reverse_complement()
revcom = revseq.reverse_complement()
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if record.seq[i:i+len(v1)] == v1:
                for j in range(10,25):
                    if record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] == v2:
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here are the first few lines of example output from the first fasta file that you linked.

seq          start    stop     gap  dir  seq
NC_003030.1  97512    97543    20   F    TTTACAAAAGAAAAGATATGCATAATGTCTTT
NC_003030.1  283173   283205   21   F    TTTACACTTATGTTTATATAGTTCAACGTCTTT
NC_003030.1  303104   303125   10   R    AAAGACGCAAGCATAATGTAAA
NC_003030.1  1127691  1127725  23   F    TTTACAAGGTGTTACTGATAGAACAAAGGGTCTTT
NC_003030.1  1301958  1301988  19   F    TTTACATAAAAAGATGATAGTTTTTGTCTTT

python 3 with redundancy

If you need sequence redundancy, this version works. I used the redundancy codes available from www.reverse-complement.com You may use the redundant characters in [ R, Y, S, W, K, M, B, V, D, H, N ] or lowercase variants. I'm not sure how to handle the '-' character at the moment but am open to suggestions. Code and results below.

Again, the algorithm isn't great given all the for loops so other answers are appreciated!

from Bio.Seq import Seq

forseq = "tttaca" # Just use python strings here. No Seq() as above
revseq = "ntnntn"
filename = "NC_003030.fasta"

revcomp = {'A':'T', 'C':'G', 'G':'C', 'T':'A',
           'R':'Y', 'Y':'R', 'S':'S', 'W':'W',
           'K':'M', 'M':'K', 'B':'V', 'V':'B',
           'D':'H', 'H':'D', 'N':'N'}
redund = {'A':['A'], 'C':['C'], 'G':['G'], 'T':['T'],
          'R':['A','G'], 'Y':['C','T'],
          'S':['C','G'], 'W':['A','T'],
          'K':['G','T'], 'M':['A','C'],
          'B':['C','G','T'], 'V':['A','C','G'],
          'D':['A','G','T'], 'H':['A','C','T'],
          'N':['A','C','G','T']}

def revcom(seq):
    return "".join(revcomp.get(base, base) for base in reversed(seq))

def seqmatch(query, reference):
    """Checks if the query sequence (perhaps a promoter with redundant bases)
    matches the reference. Both query and reference must be the same length.

    Returns true if there is a match (redundancy allowed), returns false
    otherwise.
    """
    assert len(query) == len(reference)
    for i in range(len(query)):
        if reference[i] in redund[query[i]]:
            pass
        else:
            return False
    return True

forseq = forseq.upper()
revseq = revseq.upper()
forcom = revcom(forseq)
revcom = revcom(revseq)
print("seq\tstart\tstop\tgap\tdir\tseq")
for record in SeqIO.parse(filename, "fasta"):
    for i in range(len(record.seq)):
        for v1, v2, direc in [ (forseq, revseq, "F"), (revcom, forcom, "R") ]:
            if seqmatch( v1, record.seq[i:i+len(v1)] ):
                for j in range(10,25):
                    if seqmatch( v2, record.seq[i+len(v1)+j:i+len(v2)+j+len(v2)] ):
                        print("{}\t{}\t{}\t{}\t{}\t{}".format(
                            record.id, i,
                            i + len(v1)+j+len(v2)-1,
                            j, direc,
                            record.seq[i: i+len(v1) + j + len(v2)] ))

Here is some sample output:

seq          start  stop  gap  dir  seq
NC_003030.1  733    767   23   R    GATGAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  736    767   20   R    GAAGATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  739    767   17   R    GATCAAGAAACCGATACAAACAATGTAAA
NC_003030.1  742    767   14   R    CAAGAAACCGATACAAACAATGTAAA
NC_003030.1  743    767   13   R    AAGAAACCGATACAAACAATGTAAA
NC_003030.1  2665   2697  21   R    TATTAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  2668   2697  18   R    TAAAATTTTAGAGGATGATGGTGATGTAAA
NC_003030.1  3798   3822  13   R    CACCAGAAGACTTAAAAATTGTAAA
NC_003030.1  3801   3822  10   R    CAGAAGACTTAAAAATTGTAAA


If you sort and index your BAM file:

samtools idxstats. The third column is the number of Illumina reads mapped to each PacBio read. You can then pipe that to awk: samtools idxstats foo.bam | awk '{if($1 != "*" && $3 > 0) tot += 1}END{print tot}. That will print the number of PacBio reads with at least 1 read mapped to them.
You'll want samtools depth for that

samtools depth -aa foo.bam | awk '{if($1 != lread) {if(lread != "") {print lread, 100*lcov/llen}; lread = $1; lcov=0; llen=0}; llen += 1; if($3>0) {lcov +=1}}END{print lread, 100*lcov/llen}'

That won't give you the number of illumina reads per pacbio read, but you can get that from samtools idxstats (column 3) and join the results.

BTW, your usage of "coverage" is non-standard. You want "percentage coverage" rather than "coverage", which actually means "fold-coverage" (i.e., number of reads times their average length divided by the total genome/chromosome/etc. length).

If you instead want "fold-coverage" then the easiest way is to first know the length of your reads. Presuming you have 100 base reads (for paired-end reads, this would be the same, presuming the ends don't overlap):

samtools idxstats foo.bam | awk '{print $1, $3, 100*$3/$2}'

That will give you a table of PacBio read name, number of aligned reads, fold-coverage.


There is no chromosome named chr10, but there is one named 10. The problem is that you're mixing chromosome name systems with a tool that, like most, isn't able to handle that. Remove all instances of chr in your BED file and fix things like chrM -> MT.


How to uninstall hmmpress from my conda environment?

This time conda uninstall hmmer in the root environment did the job. The sequence of commands that led me there is shown below. Thanks Devon Ryan.

conda list | grep hmmer
hmmer                     3.1b2                         3    bioconda

conda uninstall hmmer
Solving environment: done

## Package Plan ##

  environment location: /Users/user/anaconda

  removed specs:
    - hmmer

The following packages will be downloaded:

    package                    |            build
    ---------------------------|-----------------
    conda-4.5.9                |           py36_0         626 KB  conda-forge

The following packages will be REMOVED:

    barrnap: 0.9-0        bioconda
    busco:   3.0.2-py36_6 bioconda
    hmmer:   3.1b2-3      bioconda
    prokka:  1.12-3       bioconda

The following packages will be UPDATED:

    conda:   4.5.8-py36_1 conda-forge --> 4.5.9-py36_0      conda-forge
    openssl: 1.0.2o-0     conda-forge --> 1.0.2o-h470a237_1 conda-forge

Proceed ([y]/n)? y

Downloading and Extracting Packages
conda-4.5.9          |  626 KB | ################################################################################################################# | 100%
Preparing transaction: done
Verifying transaction: done
Executing transaction: done

Checking that hmmpress is not available anymore:

hmmpress
zsh: command not found: hmmpress


Bacterial genome annotation of a clinical isolate strain?

For a quick (but reliable) analysis, I'd recommend using Kallisto or Salmon to quantify isoform read counts using the transcriptome of the lab strain.

If you have a concern about transcripts that are in your sample but not in the lab strain, you can do two Trinity assemblies: one fully de-novo, and one genome-guided. These transcripts can then be used as a reference assembly for Kallisto/Salmon, and counts compared to the lab-strain results. Note that the Trinity assemblies will only apply to the particular environments that are in your samples; if a transcript isn't expressed, it won't be in the assembly.

Any transcripts in the de-novo assembly that are not in the genome-guided assembly are potentially novel to your strain. However, care must be taken in interpreting this: if there is any other contamination in the sample, those contaminant transcripts will also be assembled.

I'm not too familiar with annotation pipelines, but based on a couple of Twitter posts it looks like Prokka might be a reasonable start. Torsten Seeman (lead developer of Prokka) has done a post about alternative annotation pipelines:

http://thegenomefactory.blogspot.co.nz/2013/03/bacterial-genome-annotation-systems.html


How to convert ChiP-seq track in .bam file format into an array at fixed bin resolution?

Using deepTools to get a bin size of 10 and 20 threads:

bamCoverage -of bedgraph -bs 10 -b input.bam -p 20

You'll probably prefer the default bigWig format if you need random access or to use something like IGV.

Regarding the options to igvtools, -w is equivalent to -bs/--binSize in deepTools and denotes the resolution of the output data. -e (it's the same option in bamCoverage) denotes whether and how much to extend a read in its 3' direction. This done such that the output represents the sequences fragments. -z is specific to things like bigWig files (and apparently tdf files), where there are precomputed values that IGV will load as you zoom out (the various resolutions are "zoom levels", which is parlance taken from the bigWig format).


You can run a Molecular Dynamics simulation to see if the pose is stable and whether the protein changes/adapts it’s conformation.

You can calculate the ligand efficiency, the docking score divided by the number if heavy atoms, and rank the compounds accordingly.

You can also compare the scores of known ligands to the experimental values and see if the docking the scores are predicting known ligands or correlate with the experimental ones.

At some point you have to make a selection and purchase or synthesize the compounds. How many compounds can you test experimentally? How many can you purchase? What is the purpose of these ligands going to be? Drugs? Did you screen the correct set of ligands e.g. a library of lead compounds that are ready to purchase for cheap, that would be a good start.


cell type specific genes - heatmap using rank based approach

The figure caption says "the top 40 genes per celltype" - so we can simply use the code you provided for the top 40 genes for a single celltype, repeat this for all celltypes, and then remove any duplicate genes. Here's an example:

> expr <- matrix(runif(9*24000), nrow=24000, ncol=9)
> top_genes_each_celltype <- sapply(1:ncol(expr), function(i) order(expr[,i],decreasing =T)[1:40])
> dim(top_genes_each_celltype)
[1] 40 9
> top_genes_all_celltypes <- unique(as.vector(top_genes_each_celltype))
> length(top_genes_all_celltypes)
[1] 358

Then you can build a biclustered heatmap with expr[top_genes_all_celltypes,], for example, using heatmap2.


What i5 index should I use on the Illumina sample sheet for an unindexed p5 primer?

Illumina machines read the two indices in separate index sequencing reactions where they have sequencing barcodes specific for indices. If you have no index, then you should either not bother with creating either index read or mask that during demultiplexing. If you're mixing indexed samples with an unindexed sample, you can't really rely on any particular sequence being returned by the machine (i.e., just use the "undetermined indices" fastq file). The unindexed adapters tend to not have the full index sequencing primer binding site on them, so reads without index generally don't return much signal. In practice, they're likely to adopt barcode sequences of neighboring clusters. We see that a lot when we spike in PhiX. It has no barcode but you tend to get a small amount assigned to each sample either due to index hopping or (more likely) bleed through of signal from neighboring clusters (or small amounts of those clusters walking between nanowells and then their minority signal dominating during index sequencing).


How to extract unmatched reads using bwa and samtools?

Yes, correct but overly many steps. There is no need to convert between bam and sam. samtools can read from stdin and handles both sam and bam and samtools fastq can interpret flags, therefore one can shorten this to:
bwa mem (...options) | samtools view -o out.bam
samtools fastq -f 4 out.bam > unmatched.fastq


How to maximize fastq parsing with FastqGeneralIterator (Bio.SeqIO.QualityIO)

I had not come across FastqGeneralIterator before, but I will start using it for some of my work!

One answer is to replicate the code in the FastqPhredIterator, which does exactly what you are looking to do.

However, My personal preference would be to use the biopython infrastructure as much as possible, especially when it comes to the tricky aspects of record parsing.  If the records you want extra info on are uncommon, then using FastqPhredIterator for phred quality extraction might work for you:

import io
from Bio.SeqIO.QualityIO import FastqGeneralIterator, FastqPhredIterator

with open("input.fastq", "rU") as handle:
    for title, sequence, quality in FastqGeneralIterator(handle):
        # do something with title, sequence, and quality
        # Set special_case to True when you want more info on the record
        special_case = False
        if special_case:
            record_stream = io.StringIO("\n".join([title, sequence, "+", quality]))
            record = next(FastqPhredIterator(record_stream))
            # do something special with record

What we are doing here is combining the record strings into one and dumping them into a StringIO object, which implements the file handle duck-type. We then hand the in-memory handle to FastqPhredIterator to parse the FASTQ record again. Make sure to set special_case to True only when you want extra info on a record, as the special_case branch makes at least two extra copies of the FASTQ record.


What does it mean when PDB lists multiple organisms for the same structure?

Since this is a 40S ribosomal protein, and therefore eukaryotic-specific (as the description says) the listing of E. coli is not intuitive as you note. If you go down to "Macromolecules" and click on "Nucleic Acids/Hybrids" you can see the listed E. coli sequence is just 28bp, and all Us (poly(U) tail?). So, you can safely ignore the bacterial annotation I think.


Convert R RNA-seq data object to a Python object

A simple solution for converting Seurat objects into AnnData, as described in this vignette:

library(reticulate)
seuratobject_ad <- Convert(from=seuratobject, to="anndata", filename="seuratobject.h5ad")

Alternatively, one can use Loom, "a file format designed to store and work with single-cell RNA-seq data".

In R, install LoomR:

devtools::install_github(repo = "mojaveazure/loomR")

Convert from Seurat object to loom:

pfile <- Convert(from = pbmc_small, to = "loom", filename = "pbmc_small.loom")
pfile$close() # close loom objects when finished using them.

Then import loom object in Python using loompy, or directly as AnnData:

scanpy.api.read_loom

Alternatively, see feather.

Or export as some text format (csv, json) then import into Python.


Plotting gene models - alternatives to ggbio

It seems like the ggbio object already contains the corresponding ggplot object.
You should be able to extract the ggplot from the ggbio_obj@ggplot slot and extract your gtables.
Small example:
library(ggbio)
library(ggplot2)
library(patchwork)

fl.bam <- system.file("extdata", "wg-brca1.sorted.bam", package = "biovizBase")
wh <- as(c("chr17:41239394-41319151:+"), "GRanges")
p_cov <- autoplot(fl.bam, which = wh)
#> reading in as Bamfile
#> Parsing raw coverage...
#> Read GAlignments from BamFile...
#> extracting information...

df <-  data.frame(x= sample(1:100,size = 100,replace = TRUE))
p_hist <- ggplot(df) +
    geom_histogram(aes(x))

p_hist / p_cov
#> Error: Can't add `e2` to a ggplot object.
#> Backtrace:
#>     █
#>  1. └─patchwork:::`/.ggplot`(p_hist, p_cov)
#>  2.   └─GGally:::`+.gg`(e1, e2)
#>  3.     └─e1 %+% e2
#>  4.       └─ggplot2:::add_ggplot(e1, e2, e2name)
#>  5.         ├─ggplot2::ggplot_add(object, p, objectname)
#>  6.         └─ggplot2:::ggplot_add.default(object, p, objectname)

p_hist / p_cov@ggplot
#> `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Created on 2021-02-05 by the reprex package (v1.0.0)
Session info

sessioninfo::session_info()
#> ─ Session info ───────────────────────────────────────────────────────────────
#>  setting  value
#>  version  R version 3.6.0 (2019-04-26)
#>  os       Red Hat Enterprise Linux
#>  system   x86_64, linux-gnu
#>  ui       X11
#>  language (EN)
#>  collate  en_US.UTF-8
#>  ctype    en_US.UTF-8
#>  tz       Europe/Berlin
#>  date     2021-02-05
#>
#> ─ Packages ───────────────────────────────────────────────────────────────────
#>  package              * version  date       lib source
#>  acepack                1.4.1    2016-10-29 [1] CRAN (R 3.6.0)
#>  AnnotationDbi          1.48.0   2019-10-29 [1] Bioconductor
#>  AnnotationFilter       1.10.0   2019-10-29 [1] Bioconductor
#>  askpass                1.1      2019-01-13 [1] CRAN (R 3.6.0)
#>  assertthat             0.2.1    2019-03-21 [1] CRAN (R 3.6.0)
#>  backports              1.1.8    2020-06-17 [1] CRAN (R 3.6.0)
#>  base64enc              0.1-3    2015-07-28 [1] CRAN (R 3.6.0)
#>  Biobase                2.46.0   2019-10-29 [1] Bioconductor
#>  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor
#>  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor
#>  BiocManager            1.30.10  2019-11-16 [1] CRAN (R 3.6.0)
#>  BiocParallel           1.20.1   2019-12-21 [1] Bioconductor
#>  biomaRt                2.42.1   2020-03-26 [1] Bioconductor
#>  Biostrings             2.54.0   2019-10-29 [1] Bioconductor
#>  biovizBase             1.34.1   2019-12-04 [1] Bioconductor
#>  bit                    1.1-15.2 2020-02-10 [1] CRAN (R 3.6.0)
#>  bit64                  0.9-7    2017-05-08 [1] CRAN (R 3.6.0)
#>  bitops                 1.0-6    2013-08-17 [1] CRAN (R 3.6.0)
#>  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.0)
#>  BSgenome               1.54.0   2019-10-29 [1] Bioconductor
#>  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.0)
#>  cli                    2.0.2    2020-02-28 [1] CRAN (R 3.6.0)
#>  cluster                2.0.8    2019-04-05 [2] CRAN (R 3.6.0)
#>  colorspace             1.4-1    2019-03-18 [1] CRAN (R 3.6.0)
#>  crayon                 1.3.4    2017-09-16 [1] CRAN (R 3.6.0)
#>  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.0)
#>  data.table             1.12.8   2019-12-09 [1] CRAN (R 3.6.0)
#>  DBI                    1.1.0    2019-12-15 [1] CRAN (R 3.6.0)
#>  dbplyr                 1.4.2    2019-06-17 [1] CRAN (R 3.6.0)
#>  DelayedArray           0.12.2   2020-01-06 [1] Bioconductor
#>  dichromat              2.0-0    2013-01-24 [1] CRAN (R 3.6.0)
#>  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.0)
#>  dplyr                  0.8.5    2020-03-07 [1] CRAN (R 3.6.0)
#>  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.0)
#>  ensembldb              2.10.2   2019-11-20 [1] Bioconductor
#>  evaluate               0.14     2019-05-28 [1] CRAN (R 3.6.0)
#>  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.0)
#>  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.0)
#>  foreign                0.8-71   2018-07-20 [2] CRAN (R 3.6.0)
#>  Formula                1.2-3    2018-05-03 [1] CRAN (R 3.6.0)
#>  fs                     1.4.0    2020-03-31 [1] CRAN (R 3.6.0)
#>  GenomeInfoDb           1.22.1   2020-03-27 [1] Bioconductor
#>  GenomeInfoDbData       1.2.2    2020-04-01 [1] Bioconductor
#>  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor
#>  GenomicFeatures        1.38.2   2020-02-15 [1] Bioconductor
#>  GenomicRanges          1.38.0   2019-10-29 [1] Bioconductor
#>  GGally                 1.5.0    2020-03-25 [1] CRAN (R 3.6.0)
#>  ggbio                * 1.34.0   2019-10-29 [1] Bioconductor
#>  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.0)
#>  glue                   1.4.1    2020-05-13 [1] CRAN (R 3.6.0)
#>  graph                  1.64.0   2019-10-29 [1] Bioconductor
#>  gridExtra              2.3      2017-09-09 [1] CRAN (R 3.6.0)
#>  gtable                 0.3.0    2019-03-25 [1] CRAN (R 3.6.0)
#>  highr                  0.8      2019-03-20 [1] CRAN (R 3.6.0)
#>  Hmisc                  4.4-0    2020-03-23 [1] CRAN (R 3.6.0)
#>  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.0)
#>  htmlTable              1.13.3   2019-12-04 [1] CRAN (R 3.6.0)
#>  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.0)
#>  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.0)
#>  httr                   1.4.1    2019-08-05 [1] CRAN (R 3.6.0)
#>  IRanges                2.20.2   2020-01-13 [1] Bioconductor
#>  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.0)
#>  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.0)
#>  labeling               0.3      2014-08-23 [1] CRAN (R 3.6.0)
#>  lattice                0.20-38  2018-11-04 [2] CRAN (R 3.6.0)
#>  latticeExtra           0.6-29   2019-12-19 [1] CRAN (R 3.6.0)
#>  lazyeval               0.2.2    2019-03-15 [1] CRAN (R 3.6.0)
#>  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.0)
#>  magrittr               1.5      2014-11-22 [1] CRAN (R 3.6.0)
#>  Matrix                 1.2-17   2019-03-22 [2] CRAN (R 3.6.0)
#>  matrixStats            0.56.0   2020-03-13 [1] CRAN (R 3.6.0)
#>  memoise                1.1.0    2017-04-21 [1] CRAN (R 3.6.0)
#>  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.0)
#>  munsell                0.5.0    2018-06-12 [1] CRAN (R 3.6.0)
#>  nnet                   7.3-12   2016-02-02 [2] CRAN (R 3.6.0)
#>  openssl                1.4.1    2019-07-18 [1] CRAN (R 3.6.0)
#>  OrganismDbi            1.28.0   2019-10-29 [1] Bioconductor
#>  patchwork            * 1.0.0    2019-12-01 [1] CRAN (R 3.6.0)
#>  pillar                 1.4.3    2019-12-20 [1] CRAN (R 3.6.0)
#>  pkgconfig              2.0.3    2019-09-22 [1] CRAN (R 3.6.0)
#>  plyr                   1.8.6    2020-03-03 [1] CRAN (R 3.6.0)
#>  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.0)
#>  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.0)
#>  progress               1.2.2    2019-05-16 [1] CRAN (R 3.6.0)
#>  ProtGenerics           1.18.0   2019-10-29 [1] Bioconductor
#>  purrr                  0.3.4    2020-04-17 [1] CRAN (R 3.6.0)
#>  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.0)
#>  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.0)
#>  RBGL                   1.62.1   2019-10-30 [1] Bioconductor
#>  RColorBrewer           1.1-2    2014-12-07 [1] CRAN (R 3.6.0)
#>  Rcpp                   1.0.4.6  2020-04-09 [1] CRAN (R 3.6.0)
#>  RCurl                  1.98-1.1 2020-01-19 [1] CRAN (R 3.6.0)
#>  reprex                 1.0.0    2021-01-27 [1] CRAN (R 3.6.0)
#>  reshape                0.8.8    2018-10-23 [1] CRAN (R 3.6.0)
#>  reshape2               1.4.4    2020-04-09 [1] CRAN (R 3.6.0)
#>  rlang                  0.4.6    2020-05-02 [1] CRAN (R 3.6.0)
#>  rmarkdown              2.1      2020-01-20 [1] CRAN (R 3.6.0)
#>  rpart                  4.1-15   2019-04-12 [2] CRAN (R 3.6.0)
#>  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor
#>  RSQLite                2.2.0    2020-01-07 [1] CRAN (R 3.6.0)
#>  rstudioapi             0.11     2020-02-07 [1] CRAN (R 3.6.0)
#>  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor
#>  S4Vectors              0.24.3   2020-01-18 [1] Bioconductor
#>  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.0)
#>  sessioninfo            1.1.1    2018-11-05 [1] CRAN (R 3.6.0)
#>  stringi                1.4.6    2020-02-17 [1] CRAN (R 3.6.0)
#>  stringr                1.4.0    2019-02-10 [1] CRAN (R 3.6.0)
#>  styler                 1.3.2    2020-02-23 [1] CRAN (R 3.6.0)
#>  SummarizedExperiment   1.16.1   2019-12-19 [1] Bioconductor
#>  survival               3.1-11   2020-03-07 [1] CRAN (R 3.6.0)
#>  tibble                 3.0.0    2020-03-30 [1] CRAN (R 3.6.0)
#>  tidyselect             1.0.0    2020-01-27 [1] CRAN (R 3.6.0)
#>  VariantAnnotation      1.32.0   2019-10-29 [1] Bioconductor
#>  vctrs                  0.2.4    2020-03-10 [1] CRAN (R 3.6.0)
#>  withr                  2.4.1    2021-01-26 [1] CRAN (R 3.6.0)
#>  xfun                   0.15     2020-06-21 [1] CRAN (R 3.6.0)
#>  XML                    3.99-0.3 2020-01-20 [1] CRAN (R 3.6.0)
#>  xml2                   1.3.0    2020-04-01 [1] CRAN (R 3.6.0)
#>  XVector                0.26.0   2019-10-29 [1] Bioconductor
#>  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.0)
#>  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor


Nextflow: Can I parameterize similar processes inside a for loop?

Use channels, this is what nextflow's entire design is based on. I'm honestly surprised that spawning 5 processes with the same name in a loop even works.
Turn your parameters into regular channels, use the combine operator (https://www.nextflow.io/docs/latest/operator.html#combine) to transform the 2 parameter channels into a channel that has each combination of params/inputs. Feed this channel into the  foo  process. Afterward use branch (https://www.nextflow.io/docs/latest/operator.html#branch) to split your foo output channel into 5 different channels.
typechan = Channel.from(5,3,2,7,1)
numchan = Channel.from('prot','dna', 'rna')

cartesian = numchan.combine(typechan)

process foo {
  input:
  val numtype from cartesian
  output:
  val numtype into bar
  script:
  """
  echo $numtype
  """
}

bar.branch {
  one: it[1] == 1
  two: it[1] == 2
  three: it[1] == 3
  five: it[1] == 5
  seven : it[1] == 7
}.set {result}
result.one.view()
```


bed entries collapse based on score column

Use BEDOPS bedmap --max-element or bedmap --max-element-rand to grab the maximal-scoring element over mapped intervals, e.g.:

$ sort-bed elements.unsorted.bed5 > elements.bed5
$ bedops --merge elements.bed5 > merged.bed3
$ bedmap --max-element merged.bed3 elements.bed5 > answer.bed

The file answer.bed will contain the maximal-scored element from elements.bed5 that overlaps the merged intervals in merged.bed3.

The --max-element and --max-element-rand options differ on handling ties. The first option grabs the first match from ties by sort-bed (lexiographical) sort order; the second option grabs a random match on ties.

The option you pick would depend on what you're doing with the answer. Picking a random hit might be better if you're doing downstream statistical tests, for instance.

You can add --sci and --prec <int> options to format the score output to your desired level of precision. I think the default assumption is six digits of precision; adjust as needed.

There are also --min-element and --min-element-rand options for, conversely, reporting minimal-scoring elements.

For more information: http://bedops.readthedocs.io/en/latest/content/reference/statistics/bedmap.html#score-operations


How to estimate the phenotypic variation explained by top SNPs from a GWAS study?

Answer from @maximilian-press, converted from comment:
At a guess, the GEMMA linear mixed model is controlling for population structure and other SNPs, which your R example is not (it only fits 1 SNP?), but it's hard to be sure with this information. It would be helpful to see data samples. There is probably more context and usage information available in those resources beyond just the equations. It is especially difficult to evaluate with no terms defined as in your last example (defining the terms in the earlier ones is a good start).


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

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


You should check the User Guide:

The NCBI BLAST servers are a shared resource. We give priority to
interactive users. In order to ensure availability of the service to
the entire community, we may limit searches for some high volume
users. Interactive users of the NCBI webpages through a web browser
should not encounter problems. We will move searches of users who
submit more than 100 searches in a 24 hour period to a slower queue,
or, in extreme cases, will block the requests. To avoid problems, API
users should comply with the following guidelines:

Do not contact the server more often than once every 10 seconds.

Do not poll for any single RID more often than once a minute.

Use the URL parameter email and tool, so that the NCBI can contact you
if there is a problem.

Run scripts weekends or between 9 pm and 5 am Eastern time on weekdays
if more than 50 searches will be submitted.


A .gtf file for HG38 for HISAT2 indexing

I had the same problem and I was able to solve it with that
https://github.com/NBISweden/GAAS/blob/master/annotation/knowledge/gff_to_gtf.md


How can I locate duplicated regions in a sequence?

If you know approximately where are the duplicated regions, you can cut them out from the genome and align the two block using tools like MUMmer or Mauve.

I downloaded the sequence of ch5 you linked and extracted big regions around the genes (300bp blocks). The first block is from 70100000 - 70400000 (with SMN1 sequence at 115768 - 150839) and the second from 69200000 - 69500000 (with SMN2 sequence at 140350 - 175383).

I used progressive mauve to align the two (click on File -> align with ProgressiveMauve -> add Sequence ... -> Align).

I found that the duplicated region is almost 200k long and maps between 44690 - 225520 of the first block and 19760 - 200940 of the second block, which we still need translate to the original genome coordinates, but I believe I can leave this to you :-).

P.S. you can also somehow upload annotations to Mauve, but I was lazy, so I "hand drawn" the arrow where I knew the genes are.


Answer from @heathobrien converted from comment:
looks like that link has been dead for years:
https://www.biostars.org/p/16763/
Answer from @jason-aller converted from comment:
The biostars article mentions this:
http://baderlab.org/BINDTranslation


What is the pvalue from anova(glm1, glm2, test = “Chisq”) telling me in R when the glms have a negative binomial distribution?


bedtools unable to open file despite being tab-delimited

I suspect the reason is the "-i". The following will work,

perl -pe 's/  */\t/g' compare_me_to.txt > compare_me_to.txt2

-i will want to write onto the file and create a backup. I suspect you'll have spotted chmod +w compare_me_to.txt.

Just make sure that a double space is inserted before the * character, otherwise it will go wrong. Copying and pasting will avoid this.

This is a clearer solution because it avoids missing a space character,

perl -e 's/[ ][ ]*/\t/g' compare_me_to.txt > compare_me_to.txt2

and the solution you used .. (I did on OS X), what I'd recommend is,

perl -p -e 's/[ ][ ]*/\t/g' compare_me_to.txt # to see if it works by printing on screen without altering the file

perl -i -e 's/[ ][ ]*/\t/g' compare_me_to.txt

The  -i flag will alter the original file and store a backup without requiring renaming the file.

What's it doing?
The regex [ ][ ]* structure is saying look a space and if one or more spaces occur adjacent to each other replace everything with a single tab. [ ]* means a "space" that is present or absent and if present it can be present an infinite number of times. If the [ ]* is by itself the command risks saying replace '' (zip) or 'space' with a tab, so tabs can appear at every '' (zilch), so tabs can appear everywhere. The "*" symbol in regex work differently to the "*" wildcard character in shell.


Why do we scale the rate matrix of a substitution model to make the average rate of substitution equal 1?

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

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

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

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

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

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


Looper: can't specify the individual path for one pipeline run

Looper does not generate any of these per-sample outputs.

Looper does not have any control over where your pipeline puts its outputs.

The problem is in your pipeline, not in looper. Assuming you're using a pypiper pipeline, I suggest reading the docs; for example, the "your first pipeline" tutorial: http://pypiper.readthedocs.io/en/latest/tutorials-first.html

For real-life example, see this pipeline

outfolder = os.path.abspath(os.path.join(args.output_parent, args.sample_name))
pm = pypiper.PipelineManager(name="RRBS", outfolder=outfolder, args=args, version=__version__)


local BLAST error: BLAST Database error: Error: Not a valid version 4 database

first I recommend you to check your downloaded data with md5 to be sure there was no problem with that step. How many files did you download?
Also I was curious if you had another BLAST+ version installed on your computer. Be sure that blastn is called is a version corresponding to the database you are using.


Searching "misfolded" on the RCSB PDB turns up 64 hits. Most of these won't be relevant, but there are some interesting entries on tumor associated EGFR and misolded RNA.

In order to experimentally obtain a structure of a protein it generally has to be in a stable conformation, i.e. at low energy. Being soluble also helps. This usually corresponds to the "healthy" form, as misfolding would lead to aggregation.

That said, the structures of certain aggregates have been obtained. For example, see this recent article.

There is also the concept of mutating a protein sequence to fold to a different structure, see for instance here, though I guess you are looking for cases with the exact same sequence.


data visualization RNAseq : scaling data for PCA and cluster dendogram

Scaling (or centering) makes the genes comparable: Putting the expression levels of genes in the same scale (i.e. between 0 and 1) sustains that all of your genes contribute equally to the PCA or distance calculations. On the other hand, without this step, such calculations would be dominated by the highly/lowly expressed genes. I get "better" results with scaling.


Bootstrapping cell type identifications (scRNA-seq)

I would use permutation:
i) Get the difference that you are interested in, i.e. number/proportion of monocytes in one group vs the other.
ii) Get the relevant columns (condition, cell label, ...) from the @meta.data slot of your seurat objects and rbind().
iii) Shuffle the cell labels randomly and split the data frame into two (per condition).
iv) Compare the difference between the cell counts/proportions of interest from the two permuted data frames and save this value.
iv) Repeat iii and iv for 100s of times and calculate how many times have you observed a difference at least as extreme as your real observation (i). Once divided by the number of trials, this would give you a p-value stating how likely it is that you get the observed value.
I have to add that batch effects are severe in scRNA-seq and hence a linear model where you might be able to account for batches could be more helpful.


On Ubuntu 14.04 there is a direct route to installing RcppArmadillo. You could perhaps try the same syntax.

sudo apt-get update
sudo apt-get install r-cran-rcpparmadillo


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

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

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

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

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

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


Using PLINK to find SNPs in LD (Linkage Disequilibrium) with another set of SNPs

You can't generate LD statistics for summary data from SNPs, as the LD calculation relies on the comparison of variants at the genotype level. Plink expects data to be in a genotype matrix with either each marker for each line (probably better in your case), or each individual for each line. Plink v1.9, in particular, will accept VCF-formatted input files, when populated with 'GT' information that describes the variants for each individual.

If you're okay with using the 1000 genome population for calculating LD, plink v1.9-compatible datasets should be accessible from
here.

If not, you'll need to discuss your project with IGAP to see if you can have access to genotype-level data for LD testing. That's unlikely, because there are personal identification issues associated with non-summary data.

There's comprehensive documentation about the LD command(s) on the Plink 1.9 website here.


You want to (1) see the mean for each gene, and also to (2) calculate a ratio of expression levels of two genes, then compare it between clusters.

(1) First, notice that vlnPlot() is deprecated. Use VlnPlot().

You can try using the parameter do.sort=T:

VlnPlot(object=seuset, features.plot=c("DDB_G0267412", "DDB_G0277853"), do.sort=T)

Alternatively, you can return the ggplot2 object and then plot the means etc:

p <- VlnPlot(object=seuset, features.plot="DDB_G0277853", do.return=T)
p <- p + geom_boxplot(width=0.05)
p

(2) There are a few problems with calculating the ratio of gene expression levels. Firstly, what do you mean by gene expression level and how do you measure it? Surely you can have the read counts, but how do you interpret them? Is 30 counts a high level? Is it much more than 60 counts, or is it roughly the same? The same applies to the calculated ratios and the differences between them, even if we ignore amplification, gene length and other biases. And it is very hard to interpret ratios if the reference can also change. The problem is somewhat similar to RT-qPCR, where people use a set of reference genes whose expression has previously been shown to be invariant under the conditions.

In theory, you could use the raw counts (object@raw.data), the log + normalized counts (object@data), or the scaled counts (object@scale.data). The raw counts are biased by sequencing-depth, and the ratio of log or scaled values are not easily interpretable or intuitive. Therefore the library-size normalized (non-log) values seem to be the best.

Additionally, you could calculate the ratio of two genes either (a) for each cell (paired), or (b) for each group.
(a) is problematic, because of the zero values: you will have many NaN and Inf values, which cannot be removed without biasing the data.

Here are some functions for retrieving and plotting data from the object:

ident1 <- WhichCells(seuset, ident=1)
ident2 <- WhichCells(seuset, ident=2)

FetchData(seuset, vars.all='DDB_G0267412', cells.use=ident1, use.raw=T)

datatoplot1 <- seuset@raw.data['DDB_G0267412', ident1] / seuset@raw.data['DDB_G0277853', ident1]
datatoplot2 <- seuset@raw.data['DDB_G0267412', ident2] / seuset@raw.data['DDB_G0277853', ident2]

boxplot(datatoplot1, datatoplot2)

stripchart(datatoplot1,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=1)

stripchart(datatoplot2,
            vertical = TRUE, method = "jitter",
            pch = 21, col = "dark green",
            add = TRUE, at=2)


"DNA data" comes in several forms dependent on the technology used to produce them.
Companies like 23andMe are using SNP chips and those are available for an only rather a limited number of species. To be completely honest, I don't even know what would be the raw data they would provide, but I would suspect it could be a vcf file.
If the chosen technology would be DNA sequencing. DNA sequencing is provided by a few dozen of sequencing centres around the world, some of them associated with universities. DNA sequencing generates files that contain individual sequencing reads with assigned quality scores ("how much the sequencer trusts each base"), the typical formats are fastq (Illumina) or .bam (long read technologies). Usually, you would also need to find a suitable reference in one of the public databases, for instance, the NCBI portal for human reference. Then map your reads and do variant calling. Of course, then it depends on what would you be interested in.
If you don't care for your own DNA and you are more interested in life around you, you can check ncbi for all sorts of genomic data. If you would be interested to understand how the raw sequencing data work, ebi is a good place to look. Also, in every single paper generating/using sequencing data you there should be mention where in public databases you can find them.
Finally, let me warn you. Some of the genomic analyses are quite computationally heavy (not all though). If you don't have a specific question in mind, it would be wise to start with bacterial sequencing.
Hope you will have fun with DNA data!


Finding a single open reading frame with ribosomal binding site, using Biopython

You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.

There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading frames. You seem to be starting the reading frame at any non-stop codon instead of starting from AUGs. This may well make biological sense in some contexts (there are proteins that don't start with M) but is probably beyond the scope of what sounds like a simple assignment.


Correct for gene length or read counts in GO enrichment analysis

My understanding of the subject is that the bias of the gene length (and other bias) should be taken care when analyzing the expression and before enrichment analysis. The enrichment analysis should be done once the corrections has been performed. Because as the abstract of the GOseq paper states:

  GO analysis is widely used ... but standard methods give biased results on RNA-seq data due to over-detection of differential expression for long and highly expressed transcripts.

So first, take care of the differential expression bias by length and then use the GO to reduce complexity. How you take care of the bias in the RNA-seq data is another question. But the cqn package of Bioconductor can correct the expression by gene length and GC content. However, this correction might hurt the differential tool used (see this discussion in Bioconductor), so it might be better in some cases to use GOSeq.

Now, to the question itself:

  So I'm know unsure what to do, should I use the analysis corrected for length or read count. Or perhaps take only terms significant in both? Or only in one?

Use whatever correction method that yields better differentially expressed genes (DEG. If you find that the correction for length improves the accuracy of the predictions of DEG better than correcting by length and GC, then use that one.

Another option to obtain accurate GO terms, then you could use other testing procedures which don't rely uniquely in the Fisher test, such as those that take into account the structure of the GO graph. TopGO use this approach (note that it is a bit difficult to work with this package),this will reduce the role of the gene length bias (and probably other bias) in the resulting significant GO.


How to compute the Shannon entropy for a strand of DNA?

Why do you think the entropy of 0 is incorrect? It intuitively makes sense, as there is no uncertainty about the base at position 3, and thus there is no entropy.

However, what is plotted in a sequence logo isn't the entropy, but rather a measure of the "decrease in uncertainty" as the sequence is aligned. This is calculated by taking the entropy at this position if we randomly aligned sequences ($H_g(L)$), and subtracting from it the entropy of the alignment ($H_g(s)$):
$$
R_{sequence}=H_g(L) - H_s(L)
$$

The entropy at position 3 based on a random alignment is calculated by assuming there are 4 equally likely events (one for each base), and thus:
$$
\begin{align*}
H_g(L) & = -((1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2) + (1/4 \times -2)) \\
H_g(L) & = 2
\end{align*}
$$
Notably, $H_g(L)$ will always be 2 when dealing with nucleotide sequences.

So in your example, we have:
$$
\begin{align*}
R_{sequence}&=H_g(L) - H_s(L) \\
R_{sequence}&=2 - 0 \\
R_{sequence}&=2 \text{ bits of entropy}
\end{align*}
$$

Finally, to work out the height of each base at each position in the logo, we multiply the frequency of that base by the overall information gain $R_{sequence}$. In this case the frequency of the base A at position 3 is of course 1:
$$
\begin {align*}
H(b,l) &=f(b, l) \times R_{sequence} \\
H(A, 3) &= f(a, 3) \times 2 \\
&= 1 \times 2 \\
&= 2
\end{align*}
$$

Source: Sequence logos: a new way to display consensus sequences.
 (Thomas D. Schneider, R. Michael Stephens, 1990)

Note: I've excluded the correction factor $e(n)$ from the overall information-gain calculation for the sake of simplicity. Refer to the paper above for an explanation of this.


Score predefined ChIP-seq peaks with MACS2 or equivalent

You could try setting a high p-value threshold when you're calling peaks to retain the "non-significant" values. Something like:

macs2 callpeak -t <ChIP>.bam -c <Control>.bam -f BED -g hs -n test -B -p 0.7

You probably don't want to set the threshold to 1, since you'll likely want to avoid some noise.

But you can then use these "non-significant" values in combination with a bedtools map command, to map them to your predefined union of "actual" peaks.


Problems about value log2(IP/Input) less than zero in ChIP-seq?

Seconding what Devon says:

Use a proper normalization technique such as TMM from edgeR. For applying this to bigwig files you could follow my tutorial at biostars.org.

My two cents beyond that:

I would not use enrichment over input as input is hard to normalize and can be sparse depending on sequencing depth. Therefore, you might get biased enrichments that do not well reflect reality, especially if counts in IP or input (of both) are low. I personally only use input for peak calling, and for plots like this then compare IPs among each other.

So the strategy would be to first normalize your bigwigs (I guess you use these for the plots?) using the strategy as linked above, and then use these robust CPMs to make the plots without any input samples, only using the IPs. You can use MA-plots to explore whether samples are properly normalized. I made a little code suggestion here in this biostars post on how to do this starting basically from a matrix of raw counts using edgeR. This is basically (from the TMM normalization standpoint) the same as in the first link, but then doing pairwise MA-plots to see whether TMM properly centers the majority of points along y-axis being somewhat zero. Normalization is key in any experiment, therefore I always recommend to use these kinds of plots for diagnostics.


Apparently there is a bug in ADMIXTURE that prevents for converging when using at the same time the haploid mode and multithreading (flags --haploid and -j respectively). The problem is solved by avoiding using one of them.


Is there a Python program that is able to show convergent and divergent evolution on a phylogenetic tree?

Well in the absence of a definition of 'divergent evolution' I would simply use a definition of allopatry. If the divergence is congruent with geographic distance or geographic foci this provides some proof of divergent evolution.
In this context you are looking for something like the biogeographic BAyesian analysis such as BayArea or geographic add-ons to Beast
There isn't a singular Python library that will handle this but there is some functionality with Biopython here
from Bio import Phylo


Partially a shameless plug... I've written a genome browser, ASCIIGenome, entirely controlled through the command line. As such, it makes relatively easy to deal with questions like yours. For example, you could do:

ASCIIGenome bam/WW00292a.bam  # Load bam file(s)

Filter for reads matching the regex passed to -i option:

[h] for help: grep -i '72789|57688|70785|56299|27116|12086'

(Optionally, show sam records as strings)

[h] for help: print

Example result:

For command help type grep -h and for more filtering options see also awk -h.

Hope this helps and feel free to post questions and open issues on github!


Principle of TMT Tags in Multiplex Proteomics

TMT stands for tandem-mass-tag. The idea is to tag different samples with one of the tags then pool all samples together and run them through the mass spectrometer together. Then the reporter groups are split off and quantified (ideally) for each peptide individually. There are different methods with different numbers of isobaric labels. Isobaric means all the labels have the same mass therefore they behave exactly the same in the mass spectrometer at first. E.g. TMT11 means you have 11 different tags with identical masses, so you can pool and process 11 samples at a time.
The protein binding group attaches the label to the peptide. Then the tagged peptides from all samples are run through the chromatography and the spectrometer. At first, the MS instruments simply collects the tagged peptides with different m/z ratios at each point in time for a while. When enough peptides are collected reporter groups are cut off and the peptides then quantified. This cutting-off happens by increasing the internal energy of the collected peptides. Because the cleavage group has the weakest chemical bond it brakes first and releases the reporter group. Again, this step separates the reporter group from the rest. Finally, the cut-off reporter groups are quantified.
The different reporter groups don't have the same mass, otherwise you would not be able to differentiate between them in the MS instrument. But you need same masses before quantification to ensure that the same peptides are all collected together. Therefore, you need to balance out the different masses of the reporter groups. That is the job of the normalization group. So, the normalization group ensures that overall the tags have the same mass and the peptides with different tags are not treated differently in the MS instrument before quantification. Then for quantification the reporter groups (with different masses between each other) are cut off and the MS instrument can therefore detect and quantify them individually. There is a time dependency of the whole process which happens in multiple phases. The whole process makes it possible to compare the amounts of peptides coming from the different samples.
The chemistry in the chromatography column is quite complicated which makes it difficult to compare samples from different experiments. The presence of one peptide might enhance or suppress the signal of another peptide or impact the chromatographic retention time. Therefore, you should only compare within the same sample. More tags TMT11, TMT16 etc. allow you to run more samples in one run.
You can then use one (or more) channels (mass-tags), to tag a standard sample and then compare everything to this standard that is assumed to be identical over multiple experiments. This way you can compare even more samples than available tags. e.g. more than 11.


Your approach in general seems to be OK, but I have a few suggestions:

You are basing all your calculations on whether the pH is above or below 7. However the N and C termini have pKa values that determine whether these groups will be ionized (8 and 3, respectively). You should check whether the pH is above/below these values instead of 7.
Different amino acid side chains can also be ionized and should be accounted for when calculating the MW. Amino acid pKas can be found here
Code wise: you pass a variable unimod_location in your function but it's never used. If it's not being used then take it out.


The AC is the number of alternate alleles from called genotypes, you can see this in the vcf header:

##INFO=<ID=AC,Number=A,Type=Integer,Description="Total number of
alternate alleles in called genotypes

Assuming diploid, if we have 3 individuals, 0/0 , 0/1, 0/0, AC will be 1, because it detects 1 count of alternate allele.

If you look at your genotype calls, all are 0, so AC = 0


How can I improve the yield of MinION sequencing runs?

I attended a talk by Josh Quick at PoreCampAU 2017, in which he discussed some common barriers to getting both good sequencing yield and long read length. It mostly boils down to being more careful with the sample preparation. Bear in mind that the MinION will still sequence a dirty sample, it will just be at a reduced yield. Here are my notes from that talk:

The hardest thing about MinION sequencing is getting the sample in the first place
There are lots of different sample types and extraction methods
You can't get longer reads than what you put in; shit in leads to shit out
DNA is very stable when not moving, but very sensitive to lateral damage
The phenol chloroform method of DNA extraction is very good, and can be used with a phase-locked gel to make extraction easier
A simple salt + alcohol extraction might be the best method for extraction (because it involves the least amount of work on the DNA)
EDTA (e.g. as found in TE buffer, and many extraction kits) is not compatible with the rapid kit
The most consistently good Nanopore runs produced by Josh's lab were 1D ligations runs on R9.4; the best overall run was a phenol-chloroform extraction + rapid kit
John Tyson can tune himself out of a low-yield hole (via software)
Getting small numbers of short reads is very important
Suggested (and mostly untested) purification techniques: spin column (60-100kb); ethanol extraction (100-150kb), dialysis (150-250kb); low melting-point agarose plug (~1Mb, DNA extraction in situ)
The nanopore protocol input is in nanograms, but should really be stated as molarity; the kit expects about 0.2 pmol input
Picture molecules tethered to the surface of the membrane. You can then see that the flow cell density is independent of the sequence length
Tapestation, Qubit and Nanodrop are all a good idea; a DNA sample that can pass all three tests will work well: no short-length shoulders by Tapestation, sufficient DNA by Qubit, high purity by Nanodrop
RNA can interfere with sequencing; digesting RNA is recommended
Freezing DNA is a really bad idea. The ice crystals are very good at chopping DNA up into small pieces
DNA that is kept in the fridge is remarkably stable; can be kept for over two years (and probably indefinitely)

Update on yield:

We had a run in June 2018 which produced about 3 million reads from PCR-amplified mouse cDNA (read N50 of ~1kb), and that was in a situation where the it was the first sample prepared by someone, the input DNA amount was probably about 20% of what it should have been, and we suffered a dataloss incident. If that's what happened with a bad run, I have high hopes for our runs in the future.

Our next cDNA run in August 2018 produced about 7 million reads and a similar read N50 (i.e. total yield of about 7Gb). The run was interrupted by a Windows update during the first night (I had disabled automatic updates previously, but they had been re-enabled), and there have since been additional software and hardware updates since to improve sequencing yield. If we have a good flow cell and similar sample prep, I expect that the next cDNA run that we'll do should produce over 8 million reads, and possibly over 10 million.


What sequencing data metrics should I record?

I would say that the following metrics are useful:

Sample source (Extraction place, person, date)
Details about pre-HTS preparation (e.g. amplicon PCR, rRNA removal)
Reads per sample (+ per barcode / lane, if available)
Sequencing kit (e.g. TruSeq v3), useful for adapter trimming
Sequence length (e.g. 125bp)
Paired end / single end
If paired end, target fragment length in sample prep (e.g. 650bp)
Reference sequence metadata

Scientific name
Location of reference genome (preferably a stable web address)
Version of reference genome
Read mapping rate (per sample / barcode / lane), as a percentage of total reads
Mapping program + version

If available, potential contaminating organisms

Name / description
Location of contaminant sequence (preferably a stable web address)
Estimated mapping rate (bearing in mind that there may be some overlap with the target)
Mapping program + version

More generally, ENA / EBI has an excellent reference in sample checklists, which indicate the minimum acceptable metadata for attaching to sequence archive data:

https://www.ebi.ac.uk/ena/submit/checklists


What is the exact definition of "breakend"?

A structural variant is an inversion, deletion, duplication, or translocation. A "breakend" is one endpoint of a structural variant. This paper (Alkan et al.  Nature Review Genetics 2011) may help you understand: https://www.nature.com/articles/nrg2958
The Manta paper and documentation seem sparse on diagrams. DELLY, another structural variant caller, has diagrams in their publication. Their paper may be more of use: https://academic.oup.com/bioinformatics/article/28/18/i333/245403
Specific example: Say you're sequencing a chronic myeloid leukemia (CML) tumor. You will likely find reads that half map to the ABL1 gene on chr9 and half map to the BCR gene on chr22 (https://en.wikipedia.org/wiki/Philadelphia_chromosome). The genomic position of the location where chr9 turns to chr22 is called the breakpoint. The positions on chr9 and chr22 are the breakends.


Most universities have a ChemDraw site lincence, so you may want to check that for drawing a non–3D-embedded molecule. Alternatively there are many other programs. PyMOL is not really a good choice. Also there are online drawing tools everywhere, where you can download a SMILES (see below).
Chemical formula is the wrong way to go, unless you are "enumerating".
The chemical formula of glucose is C6H12O6. This is the same for all hexoses. What you want is a SMILES. This string gives the connectivity of the atoms. For D-glucose it's OC[C@H]1OC(O)[C@H](O)[C@@H](O)[C@@H]1O, the @ symbols are chiral centres.
In Cheminformatics making a list of all possible compound with n-carbons is a big deal, this is called enumeration and you likely don't have to reinvent wheel. Furthermore, there are datasets that are filtered by the Lipinski rule.

GDB datasets (if you don't care if you cannot buy it)
Zinc-DB  (if you care if you can buy it)
Enamine Real (if you care if you can buy it reasonably)


How to iterate protein sequences using amino acids?

In response to question 1, save the following as a python3 script (e.g., aln_parser.py):

#!/usr/bin/env python3

from Bio import AlignIO
import sys

aln = sys.argv[1]

algnmnt = AlignIO.read(open(aln, 'r'), 'clustal')

for col in range(0, algnmnt.get_alignment_length()):
    proteins = set(algnmnt[:,col])
    if len(proteins) == 1:
        print(f"Postion {col}: protein {''.join(proteins)}")

Run the script to query your alignment (clustal.aln).

python3 aln_parser.py clustal.aln

To save the output to file:

python3 aln_parser.py clustal.aln > results.txt

Requires python3 and Biopython (pip3 install biopython).


Oxford nanopore promethion; does guppy3 separate fast5 into pass and fail? If not what does?

Only if they've enabled the --qscore_filtering flag when using Guppy; this flag isn't enabled by default when running Guppy from the command line.
It's a good idea to use all fast5 reads for recalling, because accuracy improves over time, so some that previously failed may slide into the pass folder.
FYI, I'd be much happier if people asked their nanopore sequencing questions here. It's a publicly-accessible website which is indexed by search engines, so makes it easier to find answers, and easier to see the answers when found (i.e. no login needed).
Note that when recalling, you set the input and output directories (I think via the -i and -s flags). These should be different directories, and ideally different output directories for each recalling attempt, to avoid duplicate file entries. The original fast5 files shouldn't be modified by recalling as long as the input and output directories are different.


Bacterial DNA at the tail of transcriptome reads. What does that mean?

From time to time these situations are encountered, notably when PCR is being deployed in rare DNA samples.
In this instance the wiki states,

Ralstonia solanacearum is an aerobic non-spore-forming,
Gram-negative, plant pathogenic bacterium.".

I'd leave it. The simple explanation is there has been a bad mix up. The best explanation is the people loading the Illumina machine, designating the adaptors etc ... In this scenario it has nothing to do with bacterial genome integration, plant bacteria don't touch mice and aerobic bacteria wouldn't grow in cancer cells.
The second issue is what is a soil bacteria doing in a mouse cancer biological facility? For the staff this wouldn't be good news.
The singular biological  possibility is its a Pseudomonas  and some species are serious opportunistic pathogens of immune-weakened humans, but there are a lot of Pseudomonas species and this one doesn't ring any bells, BUT its mice so who knows? One species in particular is associated with hospital outbreaks, but a biomedical unit that isn't quite the same thing.

The pathogen in humans is called Pseudomonas aeruginosa. If hospitals get this infection it is BAD NEWS, for example burns victims are really susceptible and without intervention it kills them. If you have evidence that the bacterial equivalent in mice is the bacterial pathogen you have found that is a substantial step. Biological Safety Facilities (BSFs) DEFINITELY get infected. They hate it when this happens and can be very picky for fear of a contamination, but it certainly happens. If you have evidence of a major contaminent pathogen in the long run it would be of HUGE help to BSFs globally. No-one really bothers with the disease pathology of mice - because their mice and spend their entire lives going 'squeek' - but when they are models for cancer their status changes beyond recognition. So I'd erase my initial thoughts and suggest, what you have found is publishable and the application is BSF contamination - don't say that though say this could be the mouse equivalent of P. aeruginosa. Its certainly a publication and goodluck.


Do you know any virtual screening libraries with small, soluble compounds?

Lipinski and logP
Most datasets will be filtered by the Lipinski rule, plus in the some cases, such as Enamine Real whether they can be bought.
One criterion in the Lipinski rule filter is the logP, or octanol-water partition coefficient. This is super important, because a compound that cannot cross the membrane is likely of limited uses —as a consequence prodrugs are a commonly taken approach.
However, albeit obvious, not all datasets in Zinc are filtered by the Lipinski rule, they just need more browsing...
Natural?
However, you may have a valid reason to eschew logP.
Nature tries its best to have intermediate metabolites that have hydrophilic handles attaches, eg. CoA, phosphates and sugars. Therefore, if you want a natural metabolite this is a completely different question. Many of these are not very purchasable, but a PDB-derived dataset of natural ligands will likely be purchasable or easily extractable. Alternatively Zinc has several natural catalogues...
Too small
Glycerol is a 3 carbon molecule and has a MW of 90 Da. Many of the old drugs are small, but not that small. For example aspirin is 180 Da, which is twice the size of glycerol. Lipinski rule says not to go over 500 Da and 250-400 Da is a lovely range. Say you have a perfectly bound glycerol, with a perfect ligand efficiency of -1.5 kcal/mol/atom, that would mean it would bind at –9 kcal/mol in the best case scenario, i.e. it is too small to bind anything with nanomolar affinity and will actually bind a lot of things. Glycerol after all is commonly cocrystallised in protein...
The GDB datasets have enumerated compounds of a given size, but don't go that low.
Docking
There will be two further complications that will hinder the accuracy.

protonated and deprotonated atoms
polarisable groups

The former requires correct conformers, the latter requires a force field that has a Drude particle, e.g. in openMM.
Caveat against sorting by binding score
Also, it should be said that if you sort by binding score (kcal/mol or kJ/mol) you will get bigger compounds (hence the whole point of ligand efficiency) and compounds with high logP. The former is because they will have more contact points to make, the latter is because given that ∆∆G is ∆G_bound - ∆G_unbound, an unbound high-logP compound is highly unhappy in aqueous solution.


Can a data file in VCF format be converted into FASTA?

You can try gatk function FastaAlternateReferenceMaker

java -jar GenomeAnalysisTK.jar \
   -T FastaAlternateReferenceMaker \
   -R reference.fasta \
   -o output.fasta \
   -L input.intervals \
   -V input.vcf \
   [--snpmask mask.vcf]


These transformations are quite simple with the adam-submit script packaged with ADAM.

Transformation of BAM files

adam-submit -- transformAlignments sample.bam sample.alignments.adam

Transformation of VCF files

adam-submit -- transformVariants sample.vcf.gz sample.variants.adam


How can I find all the columns available in UniProt.ws package for R?

Do you have a small sample of the proteins you are looking for and what you have tried so far in R using Uniprot.ws ? (including output of it), it could be helpful to identify your problem and how to solve it.

Personally, I would rather consider Bioconductor as a platform for providing R packages dedicated to bioinformatics than a package to analyze data (but it's my opinion, maybe I'm wrong).


Stranded vs. unstranded library preparation protocols in RNAseq

Short answer: cDNA is double-stranded.

I'll try to explain the concepts surrounding this with a nanopore read, where the cDNA component of the read (from the mouse genome) is indicated in a light grey, and the adapter components are either coloured or dark grey:

The cDNA library from this RNA was prepared by adding a polyT adapter (binding to the polyA sequence of the RNA), and carrying out reverse transcription to the start of the transcript. A strand-switching process was then used to prime and extend along the opposite strand of the DNA, creating a double-stranded full-length cDNA template. This template was then amplified using PCR to increase cDNA concentration and add in barcoding sequences.

You may notice that this particular read has a polyT head, rather than a polyA tail, attached to the VNP primer in the forward orientation. I can tell by this feature that this read comes from the reverse-complement strand of a DNA transcript, and so I know that the actual RNA transcript sequence is the reverse-complement of this read. Both the forward and reverse-complement sequences could potentially come through in the sequencer, because the adapters for sequencing are attached to both sides of the cDNA template, but the transcript direction can still be determined by looking at the orientation of the adapter sequences. When mapping this sequence to the genome, only a single strand (and its associated set of kmers) needs to be considered.

When considering what happens after mapping to the genome, there are four possible match situations:

Adapter is in reverse orientation, genome is in reverse orientation
Adapter is in forward orientation, genome is in reverse orientation
Adapter is in reverse orientation, genome is in forward orientation
Adapter is in forward orientation, genome is in forward orientation

In the case of the library preparation above, situations 2 and 3 correspond to a transcript that is encoded on the forward strand of the genome, whereas situations 1 and 4 correspond to a reverse-encoded transcript. Different library preparation methods could have a different way of encoding transcript direction. See here for more details on how these match situations can be used to produce strand-specific BAM files.

If I were to use an unstranded protocol (e.g. generating a cDNA library through other means, and attaching adapters randomly to the forward or reverse transcript strand), then the adapter direction wouldn't give me information about which strand this sequence came from. I would need to consider both strands when mapping (or more correctly, when identifying the native direction of the originating transcript), because the adapter sequence has been decoupled from the transcript direction.

On Illumina systems, this strand direction calculation is performed via the chemistry on the flow cell, and the adapters are [usually] not sequenced by the sequencer at all. Adapters from one direction bind to complementary sequences on the flow cell and are extended to produce the first read, and adapters from the other direction bind to produce the second read. Even with a stranded protocol, adapters from both directions are bound (because the cDNA template is double-stranded), but the binding order (as first or second read) can be used to determine the native direction of the originating transcript.


Subset data frame by similar elements from rows and column between two data frames

You can take the PATIENT_ID from the first dataframe and use them to select from your second dataframe. Than transpose, so something like this:

ID_Names <- First_df$PATIENT_ID

rownames(Second_df) <- Second_df$Hugo_Symbol

Selected_df <- Second_df[, ID_Names]

Expected_df <- t(Selected_df)


You seem to have forgotten the space after the #. It should be as in here (Just for future readers know where to find the specs). Also it should be a single line (Maybe this error was introduced while pasting here or it has been a duplicate because there are 39 samples in the 3 line and they should be 37).

Either adding 2 more samples to the initial count and removing the newline of the third line:

39 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos

Or removing the third line:

37 2 1
# CD69neg CD69pos
CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos CD69neg CD69neg CD69neg CD69neg CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69neg CD69neg CD69pos CD69pos CD69pos CD69neg CD69pos CD69pos CD69pos CD69pos


Tool or script to parse annotated VCF files


Are there databases to annotate non-coding mutations?

VarSome allows you to annotate any variant, including non coding ones. While ACMG classification will not be available if the variant isn't in a gene's transcript (the ACMG criteria are only applicable to variants in transcribed regions), non-coding variants in transcripts will be fully annotated. For example, if you look up variant chr17-41197559-G-T which is a non-coding variant found in the UTR of one of the transcripts of BRCA1, you will get a lot of information (probably best to click on the link, the image has a lot of data):

For non-coding variants that don't fall in any transcript, you can still get frequency data, when available, as well as DANN pathogenicity predictions, GERP score, an overview of the genomic regions and other variants and/or srtuctural variants that can be found around the region, and any community annotations that may have been submitted by other users.

There is also an API you can use to automate the annotation, and VarSome Pro that lets you annotate an entire VCF file, but those are both paid services.

Disclaimer: I work for the company behind VarSome.


Best practices for deciding if two structural variants are actually the same variant?

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


Most of the reference genomes on NCBI are made of contigs or scaffolds (which is just a fancy name for contigs glued together by a bunch of unknown nucleotides). Hence, by download the whole genome, you also download all the contigs.
If you search for a specific contig (presumably a known locus) of a specific species, you can search NCBI for the species name + locus / scaffold / contig name and if it was sequenced it will pop up in the "genes" section of NCBI. For example, here is a search for Pan troglodites contig NC_036898.1 (I typed literally Pan troglodites NC_036898.1 in the search bar). You can even just type a species name, and click on "genes" and see loads of all the individual contigs sequenced.
--- edit ---
There is a bunch of ways to generate something we could call contig. However, usually what we mean is "continuous sequence assembled from sequencing reads" and it's generated from whole-genome sequencing reads. We sequence a genome, assemble the reads and what we get is a set of contigs as a representation of a haploid reference genome. Some of them will be short fragments of nucleotides with unknown functions, some of them can be long and contain hundreds of genes. The important part is that it's a continuous sequence.
With a bit more relaxed definition of "contig", we could also accept as contigs outputs of targeted sequencing (including old-fashioned sanger). It's a continuous stretch of nucleotides, usually bordered by a pair of primers. I am not 100% sure, but my best guess that most of the "gene" sequences in NCBI will be sanger-based (but I should disclaim I never investigated too much where all the NCBI sequences come from).


How to sort multiple FASTA files based on their content?

The general idea is:

iterate over all files

read the first line of each file, extract the segment bit
move the file based on that

In Python, this can be done in many ways. Here’s one way:
import os
import re

def find_segment(id):
    return re.search(r'segment (\d)', id).group(1)

for segment in range(8):
    os.mkdir(f'segment-{segment + 1}')

for entry in os.scandir():
    if not entry.is_file():
        continue

    with open(entry.path, 'r') as file:
        segment = find_segment(next(file))
    os.rename(entry.path, f'segment-{segment}/{entry.path}')

Note the absence of error handling — this code assumes that every operation succeeds; empty files, or files without “segment X” in their FASTA ID, will trigger an error, as will the usual IO failures (missing permissions, etc.).

But most people would probably use a shell scripting language instead of Python to solve this. For example, here’s something equivalent in Bash (untested):
mkdir segment-{1..8}

for file in *.fa *.fasta; do
    segment=$(sed 's/.*segment \([[:digit:]]\).*/\1/; q' "$file")
    if [[ -n $segment ]]; then
    mv "$file" "segment-$segment"
    done
done


Detect transcript isoform abundance for a specific gene in scRNA-seq

Transcript quantification is a difficult enough problem as it is, when you add the extra difficulty of going from the low read numbers available in scRNAseq if gets even more difficult.

Added to this, many scRNAseq techniques (including 10X, which is what the linked dataset is) only capture the 3' end of the gene, so you can distinguish transcript with different 3' end, but the data simply doesn't exist for doing much more with standard tools.

Lior Patchner does claims that there is benefit in looking at "transcript compatibility counts" rather than just gene counts, and that you can use his Kallisto tool for this, and the these 3' end tagging tools do give you some data for the rest of the transcript.

You could also isolate transcripts with differing three 3' ends and map reads to these using the something like the UMI-Tools pipeline.

Alternatively you could look at things on a junction level - you won't capture all junctions, but you may be able to look at differential junction usage for those you do capture using something like SUPPA or MISO. For this I would use UMI-Tools to get de-duplicated BAM files, split them on cell identity and then pass them into your tool of choice. I don't really how much you will capture.

The last thing to say is that this might be a space worth watching for future developments (cough).

Data from some other single cell methods, like SMART-seq2 or C1 will sequence the whole transcript and you may be able to apply standard tools such as Salmon, Kalisto or RSEM to get transcript usage for them, although as I said, transcript level quantification is noisey, single cell data is noisey, single-cell level transcript data will be very noisy.

Whatever, you are going to have to go back to the raw sequencing data.

Disclosure: I am the author of UMI-Tools.


Will blastn remove sequences from a search with low identity?

Yes, if a query sequence has no good alignment, blast will not return a result for it. Your command will indeed output only the best hit for each query sequence, but there is no reason to assume all of your sequences will have a match.

Those that don't have a match won't be included in the output. This sounds perfectly normal unless you are blasting the same input file against a DB created from that file. Otherwise, you don't expect all of your sequences to necessarily have a match.

Also note that since the default blastn settings are quite permissive (usually an e-value of 10, for example), don't make the mistake of assuming that a blast hit means anything. The details will always depend on the specific search you are running (I have discussed this a little in my answer here), but in the vast majority of cases, any hit with an e-value > 1 can safely be ignored. In most cases, you are only interested in e-values well under 1, actually. So you probably already have all sorts of spurious results in there anyway.

For more details, please edit your question and explain what the query and DB files are.


Reject reads with low quality bases from a Bam file through pysam

For can do this by accessing the basecall qualities from PileupRead.alignment. For example:

_BASE_QUAL_CUTOFF = 30
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            qual = read.alignment.query_qualities
            if all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual]):
                # do some stuff

This involves checking every read at every position. You could memoise this task:

def check_quality(read):
    try:
        return passed_quality[read.query_name]
    except KeyError:
        qual = read.alignment.query_qualities
        passed = all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])
        passed_quality[read.query_name] = passed
        return passed_quality[read.query_name]

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff

Or, using functools.lru_cache (thanks @winni2k):

from functools import lru_cache

@lru_cache(maxsize=1000)
def check_quality(read):
    qual = read.alignment.query_qualities
    return all([ord(c)-33 >= _BASE_QUAL_CUTOFF for c in qual])

_BASE_QUAL_CUTOFF = 30
passed_quality = dict()
with pysam.Samfile(in_file, 'rb') as bam_in:
    for column in bam_in.pileup():
        for read in column:
            if check_quality(read):
                # do some stuff


Use OpenBabel.
Something like babel -isdf file_in.sdf -opdb file_out.pdb.