# Question
Title: I have 23andme text files and would like to convert to SAM/BAM format
I would like to convert 23andme text file to NextGen to BAM file for Yfull.com to read. It is difficult to get answers on how to convert to SAM/BAM file for the 23andme text file to be converted to SAM/BAM tools

# Answer
Having done 23andme myself I can tell you that your variant file, which contains SNP genotypes, cannot be converted to a bam file. It does not contain the same information as a bam file. It may be helpful to familiarize yourself with those filetypes and the technologies used to obtain them. A SAM/BAM file contains alignments of reads obtained by sequencing. Your standard\* 23andme analysis doesn't have reads, as it uses a SNP array for determining the alleles at predefined positions. Markers which are not included in your 23andme file cannot be added by other means because those were not genotyped by the array. Therefore also Y-DNA cannot use your 23andme file or add additional markers.

\*For the sake of completeness: I am aware that under rare circumstances 23andme can perform WGS for additional variant typing, but I do not assume this is the case for OP.
> 11 votes

# Answer
If you just want to get the markers which were not reported/genotyped by 23andMe (as you mentioned in the comment), you don´t want to convert your 23andMe data to SAM/BAM (which is most probably really difficult or impossible, as explained by @Wouter De Coster), but rather to perform *imputation* on your data file. This will give you a genotype file with genotype or allele probabilities for all the markers. While it will not give you a 100% certainty for the lacking variants, it can certainly give you more information on many of them.

NB: you can imputate the Y-chromosome variants *only* if your 23andMe file contain data from Y-chromosome. If you want to know about that option, let me know.
> 3 votes

# Answer
I think the closest approximation for a SNP chip would be the intensity values. When I asked 23andMe about obtaining these, the answer was that they don't plan on adding this feature.

Some programs use the .vcf file format instead of the 23andMe format, but I think the 23andMe format is compatible with relatively more 3rd party software than some other vendors. If I remember correctly, I think DNA.land even has an intermediate .vcf file, which you can download.

I actually think it is important that 23andMe *doesn't* impute the variants in your genotype table. However, if you are trying to dig deeper into some SNP chip data, I believe you can still get genotyping for free from Genes for Good. They provide a few different "raw data" formats (such as .vcf files, with and without imputation, as well as a "23andMe format" file) and I got the impression that they would consider providing an even more raw form of data (to test the effect of generating genotypes on your own, even thought I don't think most sites will change). However, to be clear, they also do not currently provide intensity values (or some other even more raw form of data).
> 3 votes

---
Tags: sam, file-formats
---


# Question
Title: Variant calling without matched normal sample
I have WGS .bam files for 3 patients (tumour and its matched derived model namely organoid) but I don't matched normal sample. If I call variants of each patients (tumour and its matched organoid), how I can use read counts at germline heterozygous positions estimated by GATK 3.2-2 to compensate for the absence of matched normal sample? I heard people use dbsnp VCF instead of the matched normal small variant VCF. But, I don't know start from where? which GATK function does that?

> I should mention, Calling SNV and indel in many tools returns vcf but I called copy number by varscan that did not return .vcf output so I am not sure what to do for CNV

Any helps please?

# Answer
> 6 votes
If you just want to filter out calls present in dbSNP then use:

```
java -jar GenomeAnalysisTK.jar \
   -T SelectVariants \
   -R reference.fasta \
   -V patient.vcf \
   --discordance dbSNP.vcf \
   -o patient.filtered.vcf

```

`--discordance` will produce calls not present in dbSNP.

# Answer
> 0 votes
Here is a clue by scatngs

https://github.com/cancerit/ascatNgs/wiki/Human-reference-files-from-1000-genomes-VCFs

and already calculated file is here

ftp://ftp.sanger.ac.uk/pub/cancer/support-files/CPIB/ascatNgs/Human/GRCh37/

Although installing scatngs itself is a big burden

---
Tags: sam, snp, gatk, wgs, vcftools
---


# Question
Title: how do I convert a vcf file for a GWAS study (using R package vcfR)
I have a vcf file with individual level genetic data, that I read using the R package `vcfR`

```
gen <- read.vcfR('test.vcf')

Scanning file to determine attributes.
File attributes:
  meta lines: 9
  header_line: 10
  variant count: 152
  column count: 3798
Meta line 9 read in.
All meta lines processed.
gt matrix initialized.
Character matrix gt created.
  Character matrix gt rows: 152
  Character matrix gt cols: 3798
  skip: 0
  nrows: 153
  row_num: 0
Processed variant: 153
All variants processed

```

How do I convert this file to a data frame that I can use for a GWAS or Mendelian Randomization study?

I was told that this should be a data frame where each column is a genetic variant and it take values 0, 1, or 2, depending on the number of "mutations".

I never worked with genetic data before and do not know much about genetics.

# Answer
> 3 votes
If you want to do GWAS or mendelian randomisation, you can do it with Plink (v2) which should be faster than R.

---
Tags: r, vcf, gwas, data-management
---


title: Infer missing UTR features in GFF3 file

text:
# Question
Title: Infer missing UTR features in GFF3 file
I am working on a GFF file that is missing the 5'UTR and 3'UTR information. For example:

```
  ctg123 . gene      1050  9000  .  +  .  ID=gene00001;Name=EDEN
  ctg123 . mRNA      1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
  ctg123 . exon      1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
  ctg123 . exon      3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
  ctg123 . exon      5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
  ctg123 . exon      7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
  ctg123 . CDS       1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
  ctg123 . CDS       5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1

```

Is there a way to add rows in the GFF with the 5'UTR and 3'UTR ranges?

# Answer
> 3 votes
Annoyingly, you can load such a GFF3 file into R (using the GenomicFeatures package) and access the UTRs, but they don't then get saved if you use the `export` function from `rtracklayer` (this would have been the simplest solution).

Since I don't know of anything to do this off-hand, I spent a couple minutes to write one. The usage would be `addGFF3UTRs.py input.gff output.gff`. You can modify a number of the labels to appear as you'd like. This hasn't been extensively tested, but it seems to produce the correct results in a small test I just ran.

---

The code is hosted on GitHub at the link I gave above and I also include it here:

```
#!/usr/bin/env python
import argparse

def parseAttributes(kvps):
    d = dict()
    for kvp in kvps:
        k, v = kvp.split("=")
        d[k] = v
    return d

def parseGFF3(fname, topLevel="gene", secondLevel="mRNA", CDS="CDS", exon="exon"):
    f = open(fname)
    genes = dict()
    transcripts = dict()
    g2t = dict()  # A dictionary with gene IDs associated to their transcript IDs

    for line in f:
        if line.startswith("#"):
            continue
        cols = line.strip().split("\t")
        attribs = parseAttributes(cols[8].split(";"))
        ID = attribs["ID"]

        if cols[2] == topLevel:
            genes[ID] = line
            g2t[ID] = []
        elif cols[2] == secondLevel:
            parent = attribs["Parent"]
            g2t[parent].append(ID)

            # Generate a list of [line, [exon entries], [CDS entries]]
            if ID not in transcripts:
                transcripts[ID] = [line, [], []]
            else:
                transcripts[ID][0] = line
        else:
            parent = attribs["Parent"]
            if parent not in transcripts:
                transcripts[parent] = [None, [], []]
            if cols[2] == exon:
                transcripts[parent][1].append((int(cols[3]), int(cols[4]), line))
            elif cols[2] == CDS:
                transcripts[parent][2].append((int(cols[3]), int(cols[4]), line))
    f.close()

    return genes, transcripts, g2t

def findUTRs(transcripts, fivePrime = True):
    for k, t in transcripts.items():
        # If there are no CDS entries then skip
        if len(t[2]) == 0:
            continue
        # Get the strand ("." will be treated as "+")
        strand = t[0].split("\t")[6]

        # For 3' UTR, just swap the strand
        if not fivePrime:
            strand = "+" if strand == "-" else "-"

        exons = [(s, e) for s, e, _ in t[1]]
        CDSs = [(s, e) for s, e, _ in t[2]]
        exons.sort()
        CDSs.sort()
        UTRs = []
        if strand != "-":
            final = CDSs[0][0]
            for s, e in exons:
                if e < final:
                    UTRs.append((s,e))
                elif s < final:
                    UTRs.append((s, final - 1))
                else:
                    break
        else:
            final = CDSs[-1][-1]
            for s, e in exons:
                if e < final:
                    continue
                elif s > final:
                    UTRs.append((s, e))
                else:
                    UTRs.append((final + 1, e))
        t.append(UTRs)

def saveGFF(oname, genes, transcripts, g2t, fivePrime, threePrime):
    o = open(oname, "w")
    o.write("##gff-version 3\n")
    for geneID, geneLine in genes.items():
        o.write(geneLine)
        # Go through each transcript
        for transID in g2t[geneID]:
            t = transcripts[transID]
            if t[0] is None:
                continue
            o.write(t[0])
            cols = t[0].strip().split("\t")

            # Write the exons, then CDS, then 5'UTR, then 3'UTR
            for exon in t[1]:
                o.write(exon[2])
            for CDS in t[2]:
                o.write(CDS[2])
            for idx, UTR in enumerate(t[3]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], fivePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.fivePrimeUTR{};Parent={}\n".format(transID, idx, transID))
            for idx, UTR in enumerate(t[4]):
                o.write("{}\t{}\t{}\t{}\t{}\t.\t{}\t.\t".format(cols[0], cols[1], threePrime, UTR[0], UTR[1], cols[6]))
                o.write("ID={}.threePrimeUTR{};Parent={}\n".format(transID, idx, transID))
    o.close()

parser = argparse.ArgumentParser(description="Parse a GFF3 file lacking UTR entries and add them in. Note that this program assumes a reasonably formatted GFF3 file")
parser.add_argument("--fiveUTRname", default="five_prime_UTR", help="The label for 5' UTR entries (default: %(default)s)")
parser.add_argument("--threeUTRname", default="three_prime_UTR", help="The label for 3' UTR entries (default: %(default)s)")
parser.add_argument("--topLevelID", default="gene", help="The 'type' designating the top-level entry (default: %(default)s)")
parser.add_argument("--secondLevelID", default="mRNA", help="The 'type' designating the second-level entry, typically something like 'mRNA' or 'transcript' (default: %(default)s)")
parser.add_argument("--exonID", default="exon", help="The 'type' designating exons (default: %(default)s)")
parser.add_argument("--CDSID", default="CDS", help="The 'type' designating CDS (default: %(default)s)")
parser.add_argument("input", metavar="input.gff3", help="Input file name")
parser.add_argument("output", metavar="output.gff3", help="Output file name")
args = parser.parse_args()

genes, transcripts, g2t = parseGFF3(args.input, topLevel=args.topLevelID, secondLevel=args.secondLevelID, exon=args.exonID, CDS=args.CDSID)

# Add the UTRs
findUTRs(transcripts)
findUTRs(transcripts, fivePrime=False)

# Write the output
saveGFF(args.output, genes, transcripts, g2t, args.fiveUTRname, args.threeUTRname)

```

# Answer
> 4 votes
Among other things, the `canon-gff3` program in the AEGeAn Toolkit infers and prints missing features like UTRs.

```
$ canon-gff3 < genes.gff3 > genes-with-utrs.gff3

```

# Answer
> 1 votes
You can use the `add_utrs_to_gff.py` script provided by NCBI RefSeq for this purpose. Running the script on your example produces the following:

```
./add_utrs_to_gff.py example.gff3
ctg123  .  gene             1050  9000  .  +  .  ID=gene00001;Name=EDEN
ctg123  .  mRNA             1050  9000  .  +  .  ID=mRNA00001;Parent=gene00001;Name=EDEN.1
ctg123  .  exon             1050  1500  .  +  .  ID=exon00002;Parent=mRNA00001
ctg123  .  exon             3000  3902  .  +  .  ID=exon00003;Parent=mRNA00001
ctg123  .  exon             5000  5500  .  +  .  ID=exon00004;Parent=mRNA00001
ctg123  .  exon             7000  9000  .  +  .  ID=exon00005;Parent=mRNA00001
ctg123  .  CDS              1201  1500  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              3000  3902  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  CDS              5000  5300  .  +  0  ID=cds00001;Parent=mRNA00001;Name=edenprotein.1
ctg123  .  five_prime_UTR   1050  1200  .  +  .  ID=utr00002;Parent=mRNA00001
ctg123  .  three_prime_UTR  5301  5500  .  +  .  ID=utr00004;Parent=mRNA00001
ctg123  .  three_prime_UTR  7000  9000  .  +  .  ID=utr00005;Parent=mRNA00001

```

---
Tags: file-formats, gff3
---


# Question
Title: Hits in Primer-BLAST not found with programmatic blastn query
**Summary:**

I'm writing an automated tool for designing Sanger sequencing primers (other tools don't do *exactly* what I want) that will be incorporated in a production diagnostic laboratory for on-demand NGS variant verification. I'm stuck on a problem wherein off-target BLAST hits found by NCBI's Primer-BLAST aren't showing up in my programmatic BLAST query of the same genome assembly.

**Background:**

After each iteration of design with primer3, I'm using `blastn -task blastn` to search a local hg38 (patch 12) database for off-target hits. I have logic that does round-robin pairwise comparisons (i.e. compares every hit to every other hit) to determine if each hit pair could result in an amplicon (on same chromosome, opposite-stranded, perfect match at 3'-ends, plus-strand hit is upstream of minus-strand hit, are separated by \<= 1000 bp). If my logic finds a viable off-target amplicon, it triggers a total re-design with different design parameters. If no off-targets are found, those primers are finalized, so it's important that I find any and all off-targets.

To make sure that code for viable off-target amplicons is working correctly (that it's triggering re-design or reporting a design failure when appropriate), I'm using NCBI's browser-based Primer-BLAST to verify the absence of viable off-targets, which also uses genome build hg38 patch 12.

**The problem:**

Primer-BLAST is finding hits that I can't find with my code. Even when I look at every hit returned from the `blastn` query, I can't find hits that Primer-BLAST is finding. I would like to know if there is some parameter that I am missing, or parameters I could use in addition to those I'm already using, to produce the hits I am getting from Primer-BLAST.

I believe my blast database an my query fasta file are built correctly, because I am getting logical output from the query below.

**Programmatic BLAST Query Details:**

Query file (query.fasta):

```
>for_primer
GGATCTGCCTCCTTGCC
>rev_primer
CCCAACAGGGCCCTC

```

BLASTn Query Statement:

```
blastn -db blastn -db /path/to/hg38.fasta -evalue 100 -num_threads 8 -outfmt 6 std gaps nident -dust no -word_size 7 -query /path/to/query.fasta

```

The BLAST database used hg38.fa.gz (which I believe to be patch 12) from here. That was unzipped and processed with `makeblastdb`.

**Primer-BLAST Results:**

The Primer-BLAST results are here. To re-create, use the template DNA below, search database = "Genomes for selected organisms...", organism = *Homo sapiens*, max target size = 1000. The template DNA is:

```
 GAGGAGTCAGGCCTGATCTGAAGGCAGTGGGGTACCCTCTAGGGTGTCAGCAGGGCTGTGATGTCACCAGTCTTGGCTGCAGTGTGGATTAGGTAGGCGGAGGCTGGGAGGCAGTGGGAAGTGGCCGCAGGGACCGTGTGGTCAGTCAGGGTGGTGCCATTTCCTGTGTGCACACTCGGGGGCTCCCGCCACACCAGGGTGCAGACCACCCCGCCGCTCTCCCGCCAGGCAGCCTGCAGGGTGAAGGACTGACATGGCGGGAGGAGTGGCCTTCTCTGTTGGGGCGCCTGTGGCCCCCCAGCTCTTTGATAGGCAGAGCTGCAGGATCTGCCTCCTTGCCCCTGTGAGCCCAGCAGTGATTCTTGTGTCCTCCCCAGATCCATGGCTGCGGTGAGAGGGCTCGCGCCCATGCCCGGCAGAGAATCAGCCGGGAAGGGGTCCTCTATGCCGGCAGCGGGACCAAGAACGGATCCCTGGACCCAGCCAAGAGGGCCCAGCTGCAGAGGAGGCTGGATAAGAAGCTGAGTGAGCTCAGCAACCAGAGGACCAGCCGGAGGAAGGAGAGGGGGACGTGACCGGCCGCATCCTTGCACGCCCCGCGGAGCTCTCTCCATGGTAGCCCAGGGCGCTGGCAGACCATGCTCCGCCTCCACCAGGGCCACAGAGGAGCGGAGGGGCCTATGGGGGAGGAGCGGAGGGCCCTGTTGGGGAAGGTTGGGTTTTTGGACCCCAGGGATAAGCTTTTCCGATGTCACAATGTGGAGGAAAGCACCTGGGGGACAACAGTGCTCGTGCAGGTGGGGCTTGGGAAATGCACGTCCCTTCCCCTTACTCCCCGCCAAAACCCACATCCCAGCCTCTGGATCCTGGGGAAGGTTCCAGTCCCTGGAGAATACCCAGGGCCTCAAACTTGAAGTCACTCCTCCAATGTCTGGGACTTGCCAGCTCAGCCCGTTAGGATGAGGGTGCTGAGAGGAAACAGGAAACAAGACTGCGAATGGCG

```

The `blastn` output from the query above does not contain either of these two Primer-BLAST hits for 'rev\_primer' that comprise a potential amplicon from chromosome 11:

```
product length = 110
Features flanking this product:
   12759 bp at 5' side: neurexin-2-beta isoform X3
   735 bp at 3' side: RAS guanyl-releasing protein 2 isoform X5

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726566  AGT............  64726552

Reverse primer  1         CCCAACAGGGCCCTC  15
Template        64726457  ....G..........  64726471

```

Why aren't either of these two hits returned in my `blastn` query? Thanks in advance for helping me with this.

**Update:**

I noticed in IGV that off-target sequence above is lower-case in the hg38 reference genome. I've tried passing `-soft_masking false` to blastn, but got the same result.

I thought maybe the problem was due to my using hg38 from UCSC instead of that from NCBI, so I downloaded hg19 from this page and made the blast database from that file. Same result.

I then unmasked the hg38.fna file by converting all lowercase nt to uppercase, as recommended in the NCBI genome assemblies README.txt, created a BLAST database from that and ran my code again using the unmasked database. Same result, still am not detecting the off-target hit above.

# Answer
> 1 votes
I've found the problem, and it's way less interesting than I initially thought. Essentially, my "eValue" was way too low. Instead of `-evalue 100`, I'm now using `-evalue 5000` (Primer-BLAST actually uses 30,000).

For some primers, `-evalue 5000` gives ~20,000 results. It's beyond the scope of this question, but I'm now filtering on the query/subject Tm for each hit, perfect match at the query 3'-end, and number of query/subject mismatches. My pairwise hit comparisons are now detecting plausible off-target amplicons that Primer-BLAST is detecting.

Essentially, I was relying on NCBI's hit ranking system to determine what hits were relevant, but I needed cast a wider net with blastn, and then implement a rational ranking scheme of my own.

---
Tags: blast, primer, makeblastdb
---


# Question
Title: Extracting expression data in r from GEO and got a S4 object
I want to open the expression data GSE9838 dataset from GEO using GEOquery R package.

```
gset <- getGEO(filename="GSE9838_family.soft.gz")

```

The result of this code said it's S4 object of class GSE.

How can I get the expression data of GSE9838?

# Answer
> 1 votes
Well, if you use a SOFT formatted files:

```
#You can read your file
gds <- getGEO(filename="GSE9838_family.soft.gz")

# and check it
class(gds)
show(gds)

#then you usually would try to use  an exprs() accessor to retrieve the expression matrix:
exprs (gds[[1]])

# but in this case you will get an error revealing that it is a GSE datastructure
# eg. is composite structure having GSM and GPL object lists:
names(GSMList(gds))
names(GPLList(gds))

# now GPLList is probably not interesting in this case
# but you can look at the GSM-objects from GSMList:
first <- GSMList(gds)[[1]]
first
class(first)
Meta(first)

# at this point if you read the results you´ll notice that
#    title
#    [1] "X_intercept"

# and
#   An object of class "GEODataTable"
#   ****** Column Descriptions ******
#   [1] Column      Description
#   <0 rows> (or 0-length row.names)
#   ****** Data Table ******
#   data frame with 0 columns and 0 rows

```

So here you would probably wonder if the data you want is these objects. If you go back go GEO page for GSE9838, you can find that *"Raw data provided as supplementary file. Processed data is available on Series record"*. So one may want to check series record instead:

```
# you can query the GSE9838 like this:
data = getGEO("GSE9838")

# then you can check what you´ve got
names(data)
#[1] "GSE9838_series_matrix.txt.gz"
class(data$GSE9838_series_matrix.txt.gz)

# and check what´s inside:
exprs(data$GSE9838_series_matrix.txt.gz)
# but apparently it just shows one header row.

```

So it wasn´t there either. And if you check the size of the file, you´ll find that it is just 4.8 kb. My educated guess would be that the expression values are not in the dataset itself but probably in some of the supplementary files?

If you want to download them, you can:

```
filePaths = getGEOSuppFiles("GSE9838")
filePaths
untar("GSE9838/GSE9838_RAW.tar", exdir = "GSE9838/CEL")
list.files("GSE9838/CEL")

```

...what will allow you to see all the CEL files with the data.

---
Tags: r, bioconductor, microarray, geoquery
---


# Question
Title: What is the difference between a Bioinformatics pipeline and workflow?
I want to understand the difference between pipeline systems and workflow engines.

After reading A Review of Scalable Bioinformatics Pipelines I had a good overview of current bioinformatics pipelines. After some further research I found that there is collection of highly capable workflow engines. My question is then based on what I saw for argo. I would say I can be used as a bioinformatics pipeline as well.

So how do bioinformatics pipelines differ from workflow engines?

# Answer
> 18 votes
Great question! Note that from a prescriptive standpoint, the terms ***pipeline*** and ***workflow*** don't have any strict or precise definitions. But it's still useful to take a descriptive standpoint and discuss how the terms are commonly used in the bioinformatics community.

But before talking about pipelines and workflows, it's helpful to talk about programs and scripts. A ***program*** or ***script*** typically implements a single data analysis task (or set of related tasks). Some examples include the following.

* FastQC, a program that checks NGS reads for common quality issues
* Trimmomatic, a program for cleaning NGS reads
* salmon, a program for estimating transcript abundance from NGS reads
* a custom R script that uses DESeq2 to perform differential expression analysis

A ***pipeline*** or a ***workflow*** refers to a particular kind of program or script that is intended primarily to *combine other independent programs or scripts*. For example, I might want to write an RNA-seq workflow that executes Trimmomatic, FastQC, salmon, and the R script using a single command. This is particularly useful if I have to run the same command many times, or if the commands take a long time to run. It's very inconvenient when you have to babysit your computer and wait for step 3 to finish so that you can launch step 4!

So when does a program become a pipeline? Honestly, there are no strict rules. In some cases it's clear: the 10-line Python script I wrote to split Fasta files is definitely *NOT* a pipeline, but the 200-line Python script I wrote that does nothing but invoke 6 other bioinformatics programs definitely *IS* a pipeline. There are a lot of tools that fall in the middle: they may require running multiple steps in a certain order, or implement their own processing but also delegate processing to other tools. Usually nobody worries too much about whether it's "correct" to call a particular tool a pipeline.

Finally, a ***workflow engine*** is the software used to actually execute your pipeline/workflow. As mentioned above, general-purpose scripting languages like Bash, Python, or Perl can be used to implement workflows. But there are other languages that are designed specifically for managing workflows. Perhaps the earliest and most popular of these is GNU Make, which was originally intended to help engineers coordinate software compilation but can be used for just about any workflow. More recently there has been a proliferation of tools intended to replace GNU Make for numerous languages in a variety of contexts. The most popular in bioinformatics seems to be Snakemake, which provides a nice balance of simplicity (through shell commands), flexibility (through configuration), and power-user support (through Python scripting). Build scripts written for these tools (i.e., a `Makefile` or `Snakefile`) are often called pipelines or workflows, and the ***workflow engine*** is the software that executes the workflow.

The workflow engines you listed above (such as argo) can certainly be used to coordinate bioinformatics workflows. Honestly though, these are aimed more at the broader tech industry: they involve not just workflow execution but also hardware and infrastructure coordination, and would require a level of engineering expertise/support not commonly available in a bioinformatics setting. This could change, however, as bioinformatics becomes more of a "big data" endeavor.

As a final note, I'll mention a few more relevant technologies that I wasn't able to fit above.

* Docker: managing a consistent software environment across multiple (potentially dozens or hundreds) of computers; Singularity is Docker's less popular step-sister
* Common Workflow Language (CWL): a generic language for declaring how each step of a workflow is executed, what inputs it needs, what outputs it creates, and approximately what resources (RAM, storage, CPU threads, etc.) are required to run it; designed to write workflows that can be run on a variety of workflow engines
* Dockstore: a registry of bioinformatics workflows (heavy emphasis on genomics) that includes a Docker container and a CWL specification for each workflow
* toil: a production-grade workflow engine used primarily for bioinformatics workflows

---
Tags: software-usage, terminology, workflow-management
---


title: Significant gene set testing - limma

text:
# Question
Title: Significant gene set testing - limma
While I have used my title, it is more of a general statistics question regarding how these significance tests are carried out, these include `decideTests` and `mroast` / `roast` among others.

To give a brief overview of my experimental design, I am analysing bulk RNASeq datasets and there are two main factors in the design matrix: `patientID` and `casetimeViremia`; where the latter identifies whether a patient is viremia positive or negative at a particular time point (1-4). This leaves us with 5 different possibilities for `casetimeViremia`: 1-min, 2-plus, 3-min, 3-plus, 4-min.

The design matrix was created in R by running:

`model.matrix(~patientID + casetimeViremia, data = dgeObject$samples)`

Where `dgeObject` is an object of the function `DGEList`.

For the purpose of answering this biological question, we wish to conduct the differential expression testing by seeing whether all the previously mentioned `casetimeViremia` timepoints are separately and significantly expressed relative to each person's *first* timepoint, which is labelled as: 1-min.

Now, I understood that should I wish to conduct this test, I do not need to specify a contrast matrix in the `mroast` function and can just do this by calling the column number of interest.. But in that case, how does `mroast` *know* that I wish to do the contrast relative to the 1-min timepoint? And, just wondering, should I wish to define a contrast matrix, how would I do this given that there are two `colnames` missing from the design, in this case one patientID and the 1-min timepoint?

Thank you all for your time.

# Answer
> 2 votes
The `mroast` function has an argument to specify which contrast do you want to test, quoting from the help page:

> contrast
> contrast for which the test is required. Can be an integer specifying a column of design, or the name of a column of design, or a numeric contrast vector of length equal to the number of columns of design.

So if you have the contrast of 1-min timepoint positive vs 1-min timepoint negative you can use this.
Read the help pages in detail, specially the help pages in limma are quite extensive and detailed to answer most questions if not all about how to use the software.

---
Tags: r, rna-seq, statistics, limma
---


# Question
Title: What is the difference between a MeSH concept and a MeSH descriptor
Medical Subject Headings (MeSH) is a very comprehensive vocabulary to index the biomedical literature on Medline or PubMed (which is the major Biomedical literature database). MeSH is arranged into Descriptor \> Concept \> term hierarchy, but I still do not understand the difference between a MeSH descriptor and a MeSH concept.

For example, for Knee

* MeSH descriptor is D007717
* MeSH concept is M0012066 (Under the concepts section in the link)

Even after reading the description here and here, I am still not sure about the purpose of MeSH concept and MeSH descriptor.

Could anyone please shed some light on it?

# Answer
MeSH is organized into `Descriptor > Concept > Term hierarchy`, where all the synonymous `<term>`'s are grouped together under a `<concept>`. For example, consider the MeSH terms,

1. Blood-Brain Barrier
2. Hemato-Encephalic Barrier
3. Brain-Blood Barrier

All these terms are group under the umbrella concept of "Blood-brain barrier".
> 1 votes

---
Tags: ontology, literature, nomenclature
---


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

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

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

---
Tags: rna-seq, deseq2, limma
---


title: Adding a string to start of one column

text:
# Question
Title: Adding a string to start of one column
I have this dataframe

```
> head(mut)
  sample chr     pos ref alt
1  t_005   1  928006   G   A
2  t_005   1 1649842   G   T
3  t_005   1 2020408   G   A
4  t_005   1 2031677   T   A
5  t_005   1 2901878   G   A
6  t_005   1 3039338   G   A
>

```

How I can add `"chr"` to the `chr` column? I mean, instead of `"1"` I want `"chr1"`, etc.

# Answer
This should work:

```
mut$chr <- paste("chr", mut$chr, sep = "")

```
> 4 votes

---
Tags: r
---


# Question
Title: Plot a circos plot to show the consistency between 2 samples
I have called SNP and INDEL in two matched samples by *strelka* and extract this information from .vcf file and I have these

```
CHROMOSOME  POS REF ALT SAMPLE
1   928006  G   A   t_005
1   1649842 G   T   t_005
1   2020408 G   A   t_005
1   2031677 T   A   t_005

```

and

```
CHROMODOME  POSITION    REF ALT SAMPLE
1   14115   A   T   o_005
1   541052  T   C   o_005
1   1088123 T   G   o_005
1   1232501 A   G   o_005

```

I need to show if and how extend the mutations between these two samples are consistent by a circos plot but really I don't know how people do that

Something like this

or

# Answer
> 1 votes
A circos plot is most likely not the most appropriate solution here. What I would suggest is a confusion matrix, of which you can find an example here:

For every variant in your vcf you'll add a number in this matrix. One sample is the columns, the other is the lines. If your variant is homozygous in both, then you add in that square +1 (the cell with 5845 in the example).

A perfect concordant sample pair will have only variants on the diagonal.

Here is some python code to get such a matrix. It uses cyvcf2 and pandas, and expects as input a vcf file with both samples.

```
from argparse import ArgumentParser
from cyvcf2 import VCF
import pandas as pd

def main():
    args = get_args()
    confusion_matrix(args.vcf)

def confusion_matrix(vcff):
    """
    First level of the dict is the "first" call, second level is the "second" sample
    0: hom_ref
    1: heterozygous
    2: unknown/nocall
    3: hom_alt
    """
    zygosities = {0: {0: 0, 1: 0, 2: 0, 3: 0},
                  1: {0: 0, 1: 0, 2: 0, 3: 0},
                  2: {0: 0, 1: 0, 2: 0, 3: 0},
                  3: {0: 0, 1: 0, 2: 0, 3: 0},
                  }
    for v in VCF(vcff):
        zygosities[v.gt_types[0]][v.gt_types[1]] += 1
    zygs = [2, 0, 1, 3]
    df = pd.DataFrame(index=zygs, columns=zygs)
    for tr in zygs:
        for te in zygs:
            df.loc[tr, te] = zygosities[tr][te]
    df.columns = ['nocall', 'hom_ref', 'het', 'hom_alt']
    df.index = ['nocall', 'hom_ref', 'het', 'hom_alt']
    print(df)

def get_args():
    parser = ArgumentParser(description="Create confusion matrix of SNV calls")
    parser.add_argument("vcf", help="vcf containing two samples")
    return parser.parse_args()

if __name__ == '__main__':
    main()

```

---
Tags: r, mutations, wgs, software-usage
---


# Question
Title: Displaying soft-clipped nucleotides in samtools tview
There are nicer genomics visualization tools available, but the `samtools tview` command is almost always my go-to for a quick first look at read alignments. I just brought up the following locus in `tview`.

```
81055951  81055961  81055971  81055981                                                81055991  81056001  81056011  81056021
GTGGAGGTCGGCAGCAGAGCCTGTGGCAGCT**********************************************CTGAGGGTCGGCGTGGCCTCCTGGTGGAAGTTGCACCTGGTTTGCTT
...............................                                              ...............................................
...............................**********************************************..       ......................................
   ,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,       .....T................................
.....                  ,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
.......                      ,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,                                                                    ...............................................
.............                                                                .........................................A.....
........................                                                     ...............................................
...............................**********************************************..         ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,           ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************.A......                 ,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************..                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,                         ....................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,    ,,,,,,,,,,,,,,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,        ,,,,,
...............................**********************************************..                               ,,,,,,,,,,,,,,
...............................**********************************************..
...............................**********************************************..
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..................T............**********************************************...............................................
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
...............................**********************************************...............................................
...............................**********************************************...............................................
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
..............G................**********************************************...............................................
,,c,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
...............................CTAAGCAAGTAACTGTGGTGTTAGAACCTACCTCCCTCGCCGGGCA...............................................
...............................**********************************************...............................................
...............................**********************************************........................................A......
...............................**********************************************...........................................T...
,,,,,,,,,,,,,,,,,,,,,g,,,,,,,,,**********************************************,,
,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,**********************************************,,
                                                                             ...............................................
                                                                             ...............................................
                                                                             ...............................................
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,a,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,
                                                                             ,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,

```

At first glance there appears to be only a single read spanning this insertion, although the abrupt end of many reads at the insertion breakpoints made me suspicious. Indeed, upon closer inspection, 19 additional reads have primary alignments with soft clipping at the insertion breakpoints.

By default, soft-clipped portions of the reads are not displayed. **Is there any way to toggle this behavior?**

# Answer
> 1 votes
A program called gridss does this. It can be found here. It basically extracts the clipped bases and repeatedly realigns them to reference with the aligner of your choice with bwa as default. It will do this until there are no further realignments. The simple command-line for this tool after install is:

```
java -Xmx512M -cp gridss-VERSION-with-dependencies.jar gridss.SoftClipsToSplitReads I=your_input.bam O=your_output.bam REFERENCE_SEQUENCE=your_reference.fa

```

I hope this helps!

---
Tags: samtools, read-mapping, visualization
---


# Question
Title: How to convert ChiP-seq track in .bam file format into an array at fixed bin resolution?
Is there a out of the box to convert a .bam file to a .wig file at specific resolution/bin size?

I was thinking to use igvtools, but I am not sure how it works:

`igvtools count -z 5 -w 25 -e 250 alignments.bam alignments.cov.tdf hg19`

What is the diffrence between `-w,-z,-e`?

```
-z num  Specifies the maximum zoom level to precompute. The default
               value is 7 and is sufficient for most files. To reduce file
               size at the expense of IGV performance this value can be
               reduced.

 -e, --extFactor num

The read or feature is extended by the specified distance in bp prior to counting. This option is useful for chip-seq and rna-seq applications. The value is generally set to the average fragment length of the library minus the average read length.

-w, --windowSize num

The window size over which coverage is averaged. Defaults to 25 bp.

```

# Answer
> 1 votes
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).

---
Tags: chip-seq
---


# Question
Title: Clustering of gene co-expression network by igraph R package
I have constructed a gene co-expression network from RNA-seq data. The network has more than 10000 nodes and more than 1 million edges.

The network file as an edge list format of memory around 1gb which was created by calculating Pearson correlation of each gene pairs and gene pairs having correlation greater than 95% we're selected to create the edge list. The file is too large to run in cytoscape due to limitation of computer configuration.

In R, I am able to calculate the topological properties of the network. And to cluster the gene coexperssion network, I tried using different community detection methods and but I am unable to identify which community detection is good for obtaining clusters?.

After clustering how to rank the subclusters and save each subcluster in separate files.

# Answer
For the first part, do you mean that the file is too large to be run on your computer?

For the second, if I understood correctly, you can use `igraph` or `ggnet2` to color code the nodes based on the conditions you have and see if communities are formed based on condition. This would be a simple way to do this kind of thing.

I personally prefer using `ggnet2` as I found it more user friendly. What you can do if you have an adjacency matrix, say `adj_mat` is to create a color vector based on the conditions for each node in the network, where `length(color_vec)` = `ncol(adj_mat)` = `nrow(adj_mat)` \- since it is a square matrix. Following that, you could do something like this:

```
library(ggplot2)
library(network)

neto <- network(adj_mat, directed = FALSE, loops = FALSE, hyper = FALSE)
#assuming that it is undirected here and a simple network diagram.
#run ?network for more details in case your analysis is to be more complex

color_vec <- as.factor(color_vec)

col_pal <- c("blue", "red", "green") #assuming you only have 3 colors..

names(col_pal) <- levels(color_vec)
neto %v% "color" <- as.character(col_vec)

ggnet2(neto, node.size = .9, color = "color", palette = col_pal,
       edge.size = .05, edge.color = "white", mode = "adj")

```
> 2 votes

---
Tags: r, rna-seq, gene, networks, gene-expression
---


# Question
Title: Viral genome assembly using broad viral ngs pipeline?
I am trying to assemble RNA virus genome using Broad Viral NGS pipeline BROAD VIRAL NGS PIPELINE. I am two questions :

1) As this pipeline requires unaligned bam format as input, how do I convert fastq to bam ?

2) Some of the tools like assemble\_trinity require Trimmomatic clipDb. How can I get this clipDb ?

# Answer
Since you are already using the Broad Tools sets you can use Picard FastqToSam to make the conversion

As far a clipDb I am unfamiliar with that and a quick google search and look at the trimmomatic manual were unhelpful
> 4 votes

---
Tags: genome, assembly, covid-19, sars-cov-2
---


# Question
Title: Phyre2 vs ITasser, completely different models generated
Does anyone have experience generating pdb structures with Phyre and ITasser online tools. The results generated from each given the same amino acid sequence input are very different and I am wondering whether or not this is a usual experience. I know ITasser was the number 1 rated in the CASP trails, but should the results really be this disparate?

# Answer
> 5 votes
I'm less familiar with Phyre, but I-TASSER is a really sophisticated system that takes the results of a search using multiple threaders and plugs them into an ab initio simulation which tries to minimize the energy of the models by sampling many possible 3D conformations, which I don't think Phyre does.

https://en.wikipedia.org/wiki/I-TASSER#/media/File:I-TASSER-pipeline.jpg

Compare with a similar workflow schematic for Phyre:

http://www.nature.com/nprot/journal/v10/n6/images/nprot.2015.053-F1.jpg http://www.nature.com/nprot/journal/v10/n6/fig\_tab/nprot.2015.053\_F2.html

Structure prediction still has a long way to go, and you'll always get better results if there are close homologues available in the PDB, but given the consistent high performance of I-TASSER in CASP I would treat those results as more significant. That said, it can't hurt to consider multiple answers.

Edit: included blmoore's link to second half of Phyre protocol schematic

# Answer
> 1 votes
This may not be the whole answer to the question, but one key difference between Phyre and ITasser (I think I'm correct in saying), is that Phyre only returns regions of the model it considers to be sufficiently well modelled, so it may not be a full length model. ITasser on the other hand, always generates full size models.

As one of the other answers alludes to, it is very dependent on your protein of interest. If it isn't well represented in the databases, the initial threading stages will not yield good models, and the subsequent steps all suffer as a result.

---
Tags: proteins, protein-structure
---


# Question
Title: Pattern matching and removing from the name
```
colnames(countdata) = gsub("Aligned.sortedByCoord.out.bam_rem.bam", "", colnames(countdata))
colnames(countdata)
dput(colnames(countdata))
    c("BAF60A_KD_1_S13", "BAF60A_KD_2_S14", "BAF60A_RA_1_S17", "BAF60A_RA_2_S18",
    "BAF60A_VD3_1_S15", "BAF60A_VD3_2_S16", "BAF60C_KD_1_UI_S1",
    "BAF60C_KD_2_UI_S2", "BAF60C_KD_RA_1_S11", "BAF60C_KD_RA_2_S12",
    "BAF60C_KD_VD3_1_S9", "BAF60C_KD_VD3_2_S10", "VC_1_UI_S3", "VC_2_UI_S4",
    "VC_RA1_S5", "VC_RA_2_S6", "VC_VD3_1_S7", "VC_VD3_2_S8")

```

After removing the above pattern I want to remove this part from each of the colnames `_S[anynumber]` so that I can end up like this as an example `BAF60A_KD_1` but not sure how to make that pattern and match the same and remove or replace it with empty space. I tried as below which ends up like all the digits are removed.

```
colnames(countdata) = gsub("[[:digit:]_S]", "", colnames(countdata))

colnames(countdata)

c("BAFAKD", "BAFAKD", "BAFARA", "BAFARA", "BAFAVD", "BAFAVD",
"BAFCKDUI", "BAFCKDUI", "BAFCKDRA", "BAFCKDRA", "BAFCKDVD", "BAFCKDVD",
"VCUI", "VCUI", "VCRA", "VCRA", "VCVD", "VCVD")

```

Any suggestion would be helpful

# Answer
> 3 votes
Without *regex*, split on `"_"`, drop last item, then *paste* it back again:

```
sapply(strsplit(x, "_"), function(i) paste(head(i, -1), collapse = "_"))

```

Or with *regex* (stealing the code from `tools::file_path_sans_ext` function which is used to drop file extensions):

```
sub("([^_]+)\\_[[:alnum:]]+$", "\\1", x)

# [1] "BAF60A_KD_1"     "BAF60A_KD_2"     "BAF60A_RA_1"     "BAF60A_RA_2"     "BAF60A_VD3_1"    "BAF60A_VD3_2"
# [7] "BAF60C_KD_1_UI"  "BAF60C_KD_2_UI"  "BAF60C_KD_RA_1"  "BAF60C_KD_RA_2"  "BAF60C_KD_VD3_1" "BAF60C_KD_VD3_2"
# [13] "VC_1_UI"         "VC_2_UI"         "VC_RA1"          "VC_RA_2"         "VC_VD3_1"        "VC_VD3_2"

```

---
Tags: r
---


# Question
Title: taxon exclude list for searching local blast database using blastn
I am looking for a solution to exlude certain entries when searching a local blast nt database (with blastn), specifically the sequences from uncultured / environmental samples, ideally using their taxon ids. This should ideally replicate the results from the BLAST web interface, when checking the option "*Exclude Uncultured/environmental sample sequences*", see screenshot below.

Three approaches that I am aware of so far are:

(1) grep taxon IDs for any taxa with "uncultured" or "environmental" in their scientific name from names.dmp and fetch the GIs belonging to these taxa using the E-utils. Then use the -negative\_gilist parameter for blastn. However, GIs are discontinued for future records, so this method should be replaced by Acc.Ver identifiers.

(2) locally format the blast database after prepocessing the nt file by myself and exclude records belonging to the taxa list as obtained in (1).

(3) Postprocess the blast output and remove the relevant lines, but this approach might remove a lot of alignments, whereas (1) and (2) would always generate the desired number of output alignments.

So, I am looking for other suggestions, preferably without preprocessing the nt file. Ideally, one would like to have a -negative\_taxalist option for blastn, which would resemble what is possible on the NCBI BLAST web interface, where one can exclude taxa, see screenshot below.

Thanks!

# Answer
Hurray, the solution seems already to be implemented in the new BLAST database format (v5), which is usable with blast+ versions \>= 2.8.0.

It is possible to specify an exclusion as well as inclusion list of taxon ids, which, however, have to be on species level.

Species level taxon ids can be obtained from higher level taxon ids using the script get\_species\_taxids.sh, e.g.: `get_species_taxids.sh -n Enterobacterales`

See docs here: https://ftp.ncbi.nlm.nih.gov/blast/db/v5/blastdbv5.pdf

Edit:
For creating taxon ID lists locally from the NCBI taxonomy files, the program subtree can be used. See example therein for compiling a list of taxon IDs from environmental or unclassified taxa.
> 1 votes

---
Tags: blast
---


# Question
Title: Merging / Concatentating BLAST Databases
I was wondering if any one has had success merging BLAST databases from multiple references?

For a concrete example I have hg19 and a bunch of viral and bacterial genomes that all have been index separately but I want to do a single alignment to all of them. Is there a way I could say combine hg19, virus\_A, virus\_D and bacteria\_X into a single database without having to regenerate? Since the problem could get combinatorial quickly and making a db with hg19 is time consuming I was wondering if there was a way around it.

This is similar to another ticket but asking for BLAST instead of alignment indexes.

# Answer
There's an app for that, apparently.

> Often we need to search multiple databases together or wish to search a specific subset of sequences within an existing database. At the BLAST search level, we can provide multiple database names to the “-db” parameter, or to provide a GI file specifying the desired subset to the “-gilist” parameter. However for these types of searches, a more convenient way to conduct them is by creating virtual BLAST databases for these. Note: When combining BLAST databases, all the databases must be of the same molecule type. The following examples assume that the two databases as well as the GI file are in the current working directory.
>
> **Aggregate existing BLAST databases**
>
> To combine the two nematode nucleotide databases, named “nematode\_mrna” and “nematode\_genomic", we use the following command line:
>
> ```
> $ blastdb_aliastool -dblist "nematode_mrna nematode_genomic" -dbtype nucl \
>   -out nematode_all -title "Nematode RefSeq mRNA + Genomic
>
> ```

So, you can use the `blastdb_aliastool` to create merged databases just like you asked for.
> 2 votes

---
Tags: blast, makeblastdb
---


# Question
Title: Why doesn't the mouse GRCm38/mm10 refGene genome annotation file contain non-coding transcripts?
I am new to Bioinformatics and I am exploring the refGene.txt files from the UCSC genome annotation database for several species.

My question concerns the Dec. 2011 (GRCm38/mm10) assembly of the mouse genome. I have seen that the Human one (hg39) contains both coding and non-coding transcripts. But the Mouse (mm10) RefGene.txt only contains coding transcripts.

Why is that?

Source of the file: http://hgdownload.soe.ucsc.edu/goldenPath/mm10/database/refGene.txt.gz

# Answer
It wasn't the correct file. I downloaded the file again from UCSC, compared it with the file that I was using, and they were different. This new file contained both information from coding and non-coding transcripts.

To all the kind people that took interest in my problem, I am sorry.

Edited 30/12/2020: To add a bit of clarity; the link in the question is the correct link from which to download the correct file. The file with which I was working from the beginning was the wrong one, but I don't really know where I got it from the first time.
> 1 votes

---
Tags: sequence-annotation
---


# Question
Title: How to find the position-specific weight matrix (PSWM) for a list of genes having different sizes?
Wiki has the algorithm for finding position-specific weight matrix (PSWM) for a list of nucleotide sequences having an equal number of nucleotide. How do we find the PSSM for a list of nucleotide sequences having different sizes?

# Answer
> 3 votes
If your sequences are different lengths, this implies a gapped alignment with insertions/deletions. This is not possible with a regular PSSM, so you have to create a separate PSSM for each group of homologous sequences of the same sequence length.

The alternative is to use a different approach that allows gaps like profile HMMs (as implemented in HMMER3 or HH-Suite3).

---
Tags: sequence-alignment, proteins, sequence-analysis, computation
---


title: Executing PyMOL from a Shell script

text:
# Question
Title: Executing PyMOL from a Shell script
I'm trying to execute pymol from a shell script (and it is not working). I'm not executing the script on PyMOL but even if I do this, it doesn't work neither. However, if instead of a script the commands are manually written in PyMOL, everything goes well.

Here is the code:

```
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
' > scriptt.pml
pymol scriptt.pml

```

One of the problems, you will see, is that PyMOL does not remove the water molecules. Another issue is that the program does not exit. Do you have any idea what's going on there?

# Answer
Your shell script is fine, and works as expected. The reason it hangs is that you aren't telling PyMol to exit after running your script, so it's still running and waiting for the next command. Just add a `quit` at the end and it should work:

```
#!/bin/bash
pwd
echo 'fetch 1OYT
remove resn HOH
h_add elem O or elem N
select ligand, resn FSN
select prot, 1OYT and not ligand,
save lig.pdb, ligand
save prot.pdb, prot
quit
' > scriptt.pml
pymol scriptt.pml
obabel prot.pdb -xr -O prot.pdbqt
obabel lig.pdb -O lig.pdbqt
echo 'babel conversion is finished'

```
> 1 votes

---
Tags: shell, pymol
---


# Question
Title: residual Squared Coefficient of Variation (rCV²) vs Distance to Median (DM)
I am working with single cell RNA-seq data. I obtained the squared coefficient of variation (CV²) as a measure of gene expression variability:

**I want a metric to express gene expression noise that eliminates/takes into account the observed relationship between CV² of a gene and its mean expression** so I can compare this metric across genes/conditions:

So far I have found two metrics that try to tackle this problem:

* **Distance to median (DM)**:

> A median-based trend is fitted to the log-transformed CV² against the log-transformed mean. The DM is defined as the residual from the trend for each gene. This statistic is a measure of the relative variability of each gene, after accounting for the empirical mean-variance relationship

* **residual CV² (rCV²)**:

> Fit between the mean log2 expression and CV² using a gamma generalised linear model as parametrisation then computing the absolute deviation of the CV² values from the fitted CV² values.

Do you have any clue about which metric is better or if there exists any other similar metrics?

# Answer
> 2 votes
There has been a lot of work done on this problem already, particularly these two papers:

https://www.biorxiv.org/content/10.1101/576827v2

https://www.biorxiv.org/content/10.1101/574574v1

# Answer
> 1 votes
You should do both. Distance to median (DM) is the most powerful transformation IMO.

The important feature of the coefficient of variation which is more famously known as the index of dispersion (sd2/mean) is if:

* IoD is approximately 1, the distribution is random (Poison),
* IoD \<1 the distribution is a uniform distribution (very rare in biological systems)
* IoD \>\> 1 the distribution is skewed (clustered), such as a negative binomial distribution. This is common in biological systems.

An example of a uniform distribution are the distance between apple trees in an orchard. In general biological statistics IoD is rarely useful because many distributions are moderately over dispersed. This could be different however for transcriptomics. Thus if IoD results in a moderately skewed distribution you can comment on it, but its not worth taking further.

Distance to median appears to be a "non-parametric" standardisation and that is exactly what you are looking for here. Standardisation is a formal process of accommodating natural variation, e.g. time of day of the experiment, to extract the underlying trend. It is very common (in fact almost essential) in time series.

---
Tags: scrnaseq, statistics, gene-expression
---


# Question
Title: How to filter blast results from blast run against virulence factor database
I downloaded virulence factor database from here , to predict virulence genes from few genomes that my lab sequenced. I performed local blast on the database with default setting and found around 4000 genes. I want to filter this result. I am not sure how to proceed. There are evalue, pident (percent identity), qcovs (query coverage) etc that I can tweak and I did tweak. For example, if I set pident above 40% the qcovs automatically goes close to 100%. That gave me indication above 40% is pretty significant and biologically relevant. But I am not sure what percent I should consider for treating the query as having the same function as the subject -- in this case that function is virulence.

Need some guidance.

Edit1: I ran blastp

Edit2: I am blasting protein (cds predicted from genome by annotation tool prokka). My question here is what is the most reasonable way to filter putative virulence factor using blast.

# Answer
> 4 votes
The answer is it depends. 40% pident with ~100% qcovs (assuming length of protein is not very small) tells you that the virulence factor (VF) and your prokka predicted protein are likely homologs, they belong to the same (or closely related) protein family. For some research studies, you might annotate your prokka proteins as 'putative' VFs.

On the other hand, in a clinical context, even 90% pident might not be enough evidence to say these predictions are real VFs. A single amino acid change alters the protein and therefore potentially the function. You could investigate bioinformtiaclly how the changes impact the protein structure and domains etc. Ideally validate it experimentally, e.g. for a new putative E. coli toxin, clone the protein and express it in a commensal E. coli strain and test the toxicity of the new strain.

In practise, think about what question you're asking and who you want to convince. Try numerous thresholds and see how robust your results are to threshold changes (4,000 VFs seems like a lot for a 'few' genomes). You could also try the VirulenceFinder database. When annotating antibiotic resistance genes, if you use 80% pident and qcovs thresholds with a minimum protein length then most researchers will accept your predicted annotations as 'conservative', maybe the VF community has similar arbitrary but pragmatic standards.

---
Tags: blast, blastp
---


title: How I debug this code?

text:
# Question
Title: How I debug this code?
I am trying to run `CMSclassifier::classifyCMS`function on my data but I am getting this error

```
library(CMSclassifier)

> Rfcms <- CMSclassifier::classifyCMS(my_data,method="RF")[[3]]
Error in match.names(clabs, names(xi)) :
  names do not match previous names

```

This code classifies gene expression data

Code

But on example data code works

I have attached my data and example data here

Could somebody please help me in solving this error?

My data

Example data

# Answer
The problem was I should use ENTREZ rather than gene symbols
> 2 votes

# Answer
I see what you mean but I'm not sure that is the source of the error. I kinda suspect that the example data set you provided works fine. What I think is happening is the labels between the classifier, test and training data are different at some point and so classification ceases.

Your code (albeit I rarely use R) and approach and data look fine. In fact very professional, the data appears scaled (not shown). It looks good.

Perhaps what you refer to is the row name consistency in labelling against the example data set. For example row name, ARHGEF19, this is a gene ID rather than a NCBI code, viz.

```
efetch -db nucleotide -format gpc -id "LS482351" | xtract -insd gene gene

```

**Output**

```
LS482351.1  ARHGEF19

```

Elsewhere you use an AL590644.1 as the label, I assume this is the genomic region for PADI4

```
efetch -db nucleotide -format gpc -id "AL590644" | xtract -insd source organism

```

**Output**

```
AL590644.14 Homo sapiens

```

(There appears to be version 14 for this subgenomic locus).

So the labelling is a bit mixed and you can feasible convert between them using efetch and scripting a hash/dictionary. However, I don't think its the problem here because its a random forest classification and a label is a label.

It is possible the package you are using is buggy and some of the labels are not getting parsed properly e.g. with a XXX-1 tag or even if a XXX.1 tag. The example data set uses "nice whole" integers. I don't know about R but if regex approaches are used for label parsing things can easily go wrong (its not a good approach). Just a guess however and its easy to check.
> 0 votes

---
Tags: r, rna-seq, networks, modelling
---


# Question
Title: How can I ask snakemake to produce a dag where each node represents a rule?
I am using snakemake to create workflows. It is very convenient to visualise my DAG using `snakemake --dag target{sampleA,sampleB,sampleC}File | dot -Tpdf > dag.pdf`. The resulting pdf shows all the rules' dependencies to get to the target files. However, there is one node per rule/sample combination in the pdf image. When I have many samples, this creates an unreadable image with little added value. Is there an option to merge all the nodes from the same rule together so that I have the workflow for just one sample?

# Answer
I think you are looking for `--rulegraph` instead of `--dag`

```
  --rulegraph           Do not execute anything and print the dependency graph
                        of rules in the dot language. This will be less
                        crowded than above DAG of jobs, but also show less
                        information. Note that each rule is displayed once,
                        hence the displayed graph will be cyclic if a rule
                        appears in several steps of the workflow. Use this if
                        above option leads to a DAG that is too large.

```
> 7 votes

---
Tags: snakemake
---


title: log2fold change as heatmap

text:
# Question
Title: log2fold change as heatmap
From a paper:

> (D) Expression analysis of multiple lineage-specific differentiation markers in WT and PUS7-KO EBs (14 days). Heatmap shows log2 fold change (FC) PUS7-KO to WT for each individual gene (rows) in three independent experiments (columns).

They have analysed the data in EdgeR but I was wondering how did they plot fold change when they have two condition but they have three independent experiment, I haven't used EdgeR so far, use DESeq2 where I group my biological replicate. For example if I have two condition with three biological replicate then I would group them and analyse, in case of DESeq2 I would get final result with one fold change report.

How did they do in this case reporting three fold change ?

Data sets where I see three wild type and three knockout.

# Answer
> 1 votes
If I read quickly thru the paper, I see you refer to Fig 5D. In the methods they describe that they use Fluidigm (i.e. RT-qPCR) for gene expression, not RNA-seq. The data link you give is for small-RNA-seq.

---
Tags: rna-seq, figure-reproduction
---


# Question
Title: samtools view: writing to standard output failed: Broken pipe
I have the following alignment script which uses BWA:

```
#!/bin/bash
#usage: sh align_pbs.sh /scratch/waterhouse_team/benth/dbg2olc-40x/good/ output.racon2.fasta bam-soap

mkdir $3

for r1 in $(find $1 -name "*R1*.gz");
do
  base=${r1%_R1*}
  output=$(basename $(echo $r1 | sed 's/R1//g'))
  r2=$(echo $r1 | sed 's/R1/R2/g')

  #Get read group infomration:
  #Source: https://www.biostars.org/p/280837/#310132
  header=$(zcat $r1 | head -n 1)
  id=$(echo $header | head -n 1 | cut -f 1-4 -d":" | sed 's/@//' | sed 's/:/_/g')
  sm=$(echo $header | head -n 1 | grep -Eo "[ATGCN]+$")
  #echo "Read Group @RG\tID:$id\tSM:$id"_"$sm\tLB:$id"_"$sm\tPL:ILLUMINA"
  echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA"

  #cat <<EOF
  qsub <<EOF
#!/bin/bash -l

#PBS -N $output
#PBS -l walltime=48:00:00
#PBS -j oe
#PBS -l mem=70G
#PBS -l ncpus=12
#PBS -M m.lorenc@qut.edu.au
##PBS -m bea

cd \$PBS_O_WORKDIR

conda activate bwa-hic
bwa mem \
$2 $r1 $r2 \
-t 12 \
-R "$(echo "@RG\tID:${output}\tSM:${output}\tLB:${output}\tPL:ILLUMINA")" | samblaster -r | samtools view -@ 12 -h -f 0x2 -F 2316 - | samtools fixmate - - | sed 's/^/LP1-/' | samtools sort -@ 12 - -o ${3}/${output}.sorted.dedup.bam

EOF

done

```

Unfortunatelye, this pipeline cuased the following erros:

```
samblaster: Version 0.1.24
samblaster: Inputting from stdin
samblaster: Outputting to stdout
[M::bwa_idx_load_from_disk] read 0 ALT contigs
[M::process] read 1200000 sequences (120000000 bp)...
[M::process] read 1200000 sequences (120000000 bp)...
[M::mem_pestat] # candidate unique pairs for (FF, FR, RF, RR): (2, 316829, 6, 4)
[M::mem_pestat] skip orientation FF as there are not enough pairs
[M::mem_pestat] analyzing insert size distribution for orientation FR...
[M::mem_pestat] (25, 50, 75) percentile: (432, 474, 509)
[M::mem_pestat] low and high boundaries for computing mean and std.dev: (278, 663)
[M::mem_pestat] mean and std.dev: (471.59, 53.66)
[M::mem_pestat] low and high boundaries for proper pairs: (201, 740)
[M::mem_pestat] skip orientation RF as there are not enough pairs
[M::mem_pestat] skip orientation RR as there are not enough pairs
[M::mem_process_seqs] Processed 1200000 reads in 615.287 CPU sec, 58.490 real sec
samblaster: Loaded 19 header sequence entries.
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting
sed: couldn't write 567 items to stdout: Broken pipe
[E::bgzf_flush] File write failed (wrong size)
samtools fixmate: Couldn't write to output file: Broken pipe
[E::bgzf_close] File write failed
[bam_mating] error while closing output file
samtools view: writing to standard output failed: Broken pipe
samtools view: error closing standard output: -1
samblaster: Unable to write to output file.

```

Thank you in advance,

# Answer
```
[W::sam_read1] Parse error at line 1
samtools sort: truncated file. Aborting

```

SAMtools sort has been unable to parse its input, which it thought was SAM (mostly because it couldn't be recognised as another format e.g. BAM).

This is because `sed 's/^/LP1-/'` is putting `LP1-` at the front of every line. If the output of samtools fixmate is SAM, then this LP1 is garbling the SAM header lines. If the output of samtools fixmate is BAM (as it probably is), this LP1 is completely garbling the data.

What were you trying to achieve with this `sed` command?
> 3 votes

---
Tags: samtools, bwa
---


# Question
Title: Voom transformation of RNA seq raw counts data
I have RNA seq data of raw gene counts that I want to transform for linear modelling. I am trying to voom transform, to do a weighted analysis.

Data frame: prac\_count\_10

```
gene         S1     S2    S3   S4    S5
ENG0000456   0      1     10   145   24
ENG0000458   7      2     9    0     0
ENG0000657   76     12    56   10    2
ENG0000689   0      0     0    3     5

```

Code:

```
prac_prac_10 <- voom(prac_count_10[,2:5], design=NULL, lib.size=NULL,
  normalize.method = "none", span = 0.5

```

However this doesn't give an output with the gene sample names and logCPM values. Would like an output with gene id's and logCPM values for linear modelling.

Also tried:

cpm <- cpm(prac_count_10)
lcpm <- cpm(prac_count_10, log=TRUE)

````````````

# Answer
Try to use `DElist()` function before you transform, and also make rownames first.

```
counts <- prac_count_10[,-1]

rownames(counts) <- prac_count_10[,1]

DGE1 <- DGEList(counts)

```

And then continue with your voom...
> 1 votes

---
Tags: r, limma
---


# Question
Title: What is the best way to account for GC-content shift while constructing nucleotide-based phylogenetic tree?
Let's say I want to construct a phylogenetic tree based on orthologous nucleotide sequences; I do not want to use protein sequences to have a better resolution. These species have different GC-content.

If we use a straightforward approach like maximum likelihood with JC69 or any other classical nucleotide model, conserved protein coding sequences of distant species with similar GC-content will artificially cluster together. This will happen because GC-content will mainly affect wobbling codon positions, and they will look similar on the nucleotide level.

What are possible ways to overcome this? I considered the following options so far:

1. Using protein sequence. This is possible of course, but we lose a lot of information on the short distance. Not applicable to non-coding sequences.
2. Recoding. In this approach C and T can be combined into a single pyrimidine state Y (G and A could be also combined in some implementations). This sounds interesting, but, first, we also lose some information here. Mathematical properties of the resulting process are not clear. As a result, this approach is not widely used.
3. Excluding third codon position from the analysis. Losing some short-distance information again. Also, not all synonymous substitution are specific to the third codon positions, so we still expect to have some bias. Not applicable to non-coding sequence.

It should be possible in theory to have a model which allows shifts in GC-content. This will be a non time-reversible Markov process. As far as I understand there are some computational difficulties estimating likelihood for such models.

# Answer
There are models that take into account compositional heterogeneity both under the maximum likelihood and Bayesian frameworks. Although the substitution process is not time-reversible, the computations are simplified by assuming that the instantaneous rate matrix can be decomposed into an "equilibrium frequency vector" (non-homogeneous) and a symmetric, constant exchange rate matrix.

I guess all your suggestions are also valid, and I remember recoding being used successfully to reduce the GC-content bias (examples in the references above and here).
> 5 votes

# Answer
The following 2004 paper describes a way to model compositional changes across the tree, in a Bayesian framework: https://doi.org/10.1080/10635150490445779

A python package implementing this ("p4"), and improvements added along the years, is available here: https://github.com/pgfoster/p4-phylogenetics

To get started, you may find useful examples here: http://p4.nhm.ac.uk/scripts.html

This has been used in a few large-scale phylogenetic analyses.
> 3 votes

# Answer
The answer is the logDet algorithm was constructed to overcome GC% clustering.

At the time it was devised only a distance method was available/implemented, so it wasn't very powerful. The posts here imply that a Bayesian or ML approach are available and these do tightly stick to the model.

Original publication here
> 1 votes

---
Tags: phylogenetics, phylogeny, nucleotide-models
---


# Question
Title: Is there a way to assemble contigs starting from a specific sequence?
My work involves searching for marker genes/fragments in metagenomic databases (like the Sequence Read Archive). Once I find these sequences, I would like to know more about the neighboring genomic region.

Is there a way I could assemble only sequences that create a contig which contains my region of interest? Contigs which don't contain this region are not useful to me. My organism of interest might represent a minority of the metagenome, and assembling everything in the dataset would use a lot of computing power.

# Answer
> 2 votes
This is one of the primary use cases for which spacegraphcats (preprint and code) was designed. The "neighborhood queries" discussed in the paper sound particularly relevant. I don't have any personal experience with spacegraphcats, but the run guide provides some examples of how to index the complete data set and how to query sequences of interest.

---
Tags: sequence-alignment, fasta, assembly, metagenome
---


title: Viral genome finishing

text:
# Question
Title: Viral genome finishing
I have assembled poxvirus genome using Ray. The assembly is good. Out of several thousand contigs I got, I was able to get one scaffold using Contiguator tool, which is about 90% of my genome. I have well annotated reference sequence. If there any good tool for filling gaps and checking assembly quality?

Things I've tried so far:

1. IMAGE from PAGIT. It asks for ace file which I guess is obtained from Newbler. So, I am not able to use this.
2. gapfill\_py script from broad: There are some gaps still missing.

I came to know that Consed is good. But, I have windows which is not compatible for Consed.

I checked my BAM file with Tablet. There are reads that cover almost entire gene. Are there any good tools or manual ways to finish my assembly so that I can get the full genome?

# Answer
Assuming that gap filling is a domain-independent and works the same for viruses as for bacteria or eukarytes.

There are two real options I am aware of: Sealer a successor of GapFiller from ABySS package and GapCloser from SOAPdenovo2 package. According to Sealer paper it's better, but I am not sure if any indepented benchmarking have confirmed their conclusion. On my insect data GapCloser worked better than GapFiller (back then I have not tried to compare it to Sealer).
> 2 votes

---
Tags: genome, assembly, covid-19, sars-cov-2
---


title: VEP output SIFT_score unclear

text:
# Question
Title: VEP output SIFT_score unclear
We have been experimenting with VEP (Variant Effect Predictor). One of the meta data attributes that we are interested in is the SIFT score, indeed when we apply the dbNSFP plug we get a column containing the scores (named SIFT\_score). However, I don't understand why there are sometimes dots or multiple values in the fields. For example, the gene ENSG00000196924 below has 5 transcripts:

The SIFT\_score column contains several values, not 1 per transcript/rs-number...

Here is another example that confuses me (I added the SIFT\_pred column this time):

There are two mutations, the lower one can be expressed in 4 transcripts and I thus understand that there can be 4 SIFT scores, but why are all for in given in every row? Is the first one the SIFT\_score for the first transcript?

One last example, again 4 transcripts, but now 2 of the scores are dots, what does that mean?

I have been looking for quite some time now how to interpret this data, any help is appreciated.

# Answer
> 3 votes
The dbNSFP plugin from VEP accesses tables of data for each variant from dbNSFP and pulls out the values. dbNSFP provide their SIFT scores in that format: a score for every transcript affected by the variant, all on one line. The lookup is just for the variant, not for the variant/transcript combo, so they provide scores for all variant/transcript combos. You can also get a column that gives you a list of the transcripts or proteins (Ensembl\_transcriptid or Ensembl\_proteinid) in order so you know which score goes with which transcript.

A better way to get SIFT scores with VEP is to get them directly from VEP, rather than using dbNSFP. This will get the SIFT score that goes with the transcript on the line with the relevant transcript.

# Answer
> 3 votes
The first gene you mention, ENSG00000196924, actually has 6 transcripts (link to the VarSome.com page of variant rs371839875), not 5. It's just that one of them is non-coding:

So the Sift scores you see are indeed one per transcript, it's just that there are 6 because dbNSFP also includes a score for the non-coding transcript of the gene.

The dots are just there as placeholders, they mean there was no value associated with that transcript. Many tools will show some sort of symbol instead of an empty field both for clarity and for practical technical reasons.

Visiting the variant's page on VarSome gives you a clearer picture since we collapse the identical scores and also include the converted rankscore provided by dbNSFP so you can have a single number for your variant:

---

**Disclaimer**: I work for the company behind VarSome, but it's a free tool. You need to pay to annotate VCF files (unlike the 100% free VEP), but it's free to use as a lookup tool for single variants.

---
Tags: ngs, variant-calling
---


# Question
Title: Is there a tool to convert a variableStep .wig to fixedStep one?
I have a .wig file in this format:

```
variableStep chrom=chr1 span=20
59965   7
59985   10
60005   10
60025   10

```

and I would like to format it to something like this format:

```
#bedGraph section chr16:0-103584
chr16   0       59901   0
chr16   59901   59927   1
chr16   59927   59931   2

```

Please assume, I am working with a reference hg19 genome, in case the info isn't deducible from the `variableStep` .wig file. I couldn't find anything in WiggleTools

# Answer
When I had to do something similar previously, I had used the following method using `sed` and `awk`. It is not as robust in terms of error detection and correction as the python solution but it is a simple one-liner.

```
$ cat temp.wig
variableStep chrom=chr1 span=20
59965   7
59985   10
60005   10
60025   10

$ sed -r 's/variableStep chrom=(.*) span=([0-9]*)/\1\t\2/g' temp.wig \
        | awk 'BEGIN{FS="\t";OFS="\t"}{if ($1~/^chr/) {c=$1 ; s=$2} else {print c,$1,$1+s,$2}}' \
    > temp.bedgraph

$ cat temp.bedgraph
chr1    59965   59985   7
chr1    59985   60005   10
chr1    60005   60025   10
chr1    60025   60045   10

```
> 1 votes

# Answer
This is a temporary solution I came up with that worked for my specific case:

```
def convert_line(line, chrom, span, start=None, fixedStep=None, count=None):
        new_line=''
        if fixedStep is None or count is None or start is None:
                line_arr  = line.split()
                start = line_arr[0]
                end   = int(line_arr[0]) + int(span)
                val   = line_arr[1]
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)
        else:
                line_arr  = line.split()
                val   = np.float(line_arr[0])
                start = int(start) + int(fixedStep)*int(count)
                end   = start + int(fixedStep)
                new_line = '{0}\t{1}\t{2}\t{3}\n'.format(chrom, start, end, val)

        return new_line

def convert_file(input_file, output_file):
        chrom=None
        span=None
        fixedStep=None
        count=None
        start=None
        with open(input_file, 'r') as fin, open(output_file, 'w') as fout:
                for line in fin:
                        if len(line)<1 or line[0] == '#':
                                continue
                        elif 'variableStep' in line or 'fixedStep' in line:
                                for key in line.split():
                                        if 'chrom=' in key:
                                                chrom = key.split('chrom=')[-1]
                                        if 'span=' in key:
                                                span = key.split('span=')[-1]
                                        if 'step=' in key:
                                                fixedStep = key.split('step=')[-1]
                                        if 'start=' in key:
                                                start = key.split('start=')[-1]
                                                count = 0
                        else:
                                line_out = convert_line(line, chrom, span, start, fixedStep, count)
                                fout.write(line_out)
                                if count is not None:
                                        count += 1

```

input:

```
variableStep chrom=chr16 span=20
59965   7
59985   10
60005   10

```

output:

```
chr16   59965   59985   7
chr16   59985   60005   10
chr16   60005   60025   10
chr16   60025   60045   10

```
> 2 votes

---
Tags: python
---


title: A GTF annotation for GRCH375d version

text:
# Question
Title: A GTF annotation for GRCH375d version
I have been given a bunch of .bam files, that come from the alignment on GRCh37d5 of human genome by Star tool. I want to extract raw read counts from these bam files by FEATURECOUNTS tool so I need a GTF. I googled a lot but I did not find GTF for GRCh37d5. I used another GTFs and I obtained weird features. Please one of you help me in finding this specific GTF. I converted a bam to fastq and I used Salmon to quantify gene expression by a transcripome.fasta but rate of mapping was 0.06

# Answer
> 2 votes
The most recent annotation of GRCh37 assembly by NCBI RefSeq is here: https://ftp.ncbi.nlm.nih.gov/genomes/Homo\_sapiens/GRCh37.p13\_interim\_annotation/

You can download the annotation in GFF3 format from that path and use it with featurecounts tool by first converting it to the SAF format using the commands shown below:

```
$ zgrep -v '^#' annotation.gff.gz \
    | awk 'BEGIN{FS="\t";OFS="\t"}($3=="exon"){print $1,$4,$5,$7,$9}' \
    | sed -r 's/ID.*GeneID:([0-9]*).*/\1/g' \
    | awk 'BEGIN{FS="\t";OFS="\t"}{print $5,$1,$2,$3,$4}'
    > annotation.saf

```

You should confirm that the chromosome sequence identifiers in the BAM file and the SAF file match. NCBI annotation has identifiers in the accession.version format whereas UCSC uses `chr#` format and Ensembl uses `#`. For this, you can use the `-A` parameter from `featureCounts`. Assuming that your BAM file uses UCSC-style identifiers, you can generate a 'chromosome name alias file' that can be used with featureCounts as follows:

```
## get NCBI assembly report for GRCh37
wget -nv ftp://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/001/405/GCF_000001405.24_GRCh37.p12/GCF_000001405.24_GRCh37.p12_assembly_report.txt

## generate chr_alias.csv file
grep -v '^#' GCF_000001405.24_GRCh37.p12_assembly_report.txt \
    | dos2unix \
    | awk 'BEGIN{FS="\t";OFS=","}($10!="na"){print $10,$7}' \
    > chr_alias.csv

```

---
Tags: rna-seq, sequence-alignment, sequence-annotation, gtf
---


# Question
Title: Are phylogenetic tree construction algorithms any different than general clustering algorithms?
Are phylogenetic tree construction algorithms any different from clustering algorhithms? I suspect the answer is no.

Of course phylogenetic tree construction uses biological knowledge, e.g special distance metrics, but does it brings anything new to clustering algorithm (e.g hierarchical, neighbor-joining etc).

# Answer
The goal of a phylogeny is to estimate the "expected" number of mutations between all taxa in the analysis and their hypothetical common ancestors. A cluster-analysis will only identify the "observed" mutations and "expected" and "observed" mutations can be majorly different due to the major artefact of reversion mutation. This is particularly true of nucleotide phylogenies.

The key difference between clustering algorithms based on a "uncorrected" distance matrix and phylogeny is the latter is based on an explicit model to accommodate reversion mutations. The real problem is that there are only 4 bases, so randomly there is 1/4 chance that a mutation at a given position will mutate back to the original, e.g. A-\>C-\>A. Observed mutational differences = 0, expected (real) mutational differences = 2. What phylogeny is concerned about is reconstructing that "2". The problem is significant because almost every gene has regions of rapid mutations and regions of low mutations.

The primary way of doing this is via the Jukes-Cantor correction and is foundational in all nucleotide trees expect "p-distances". If nucleotide divergence is lower than 75% then it is possible to estimate the expected number of mutations using the observed via the JC correction. Moreover, when combined with a method of estimating rate variation within a gene, (usually the discrete gamma distribution), JC correction is very effective at recovering the "true tree". This is because fast evolving homologous sites are grouped together - big correction via JC, slow evolving sites are grouped together - small correction via JC. Other approaches to improving the JC correction are via identifying the bias between purine to purine mutations and pyrimidine to pyrimidine and purine to pyrimidine mutations.

The importance of JC was demonstrated by simulation studies of 4 taxa ((a,b), (c,d)), if two strains evolved very quickly when there sister lineages evolve slowly a clustering algorithm will report the fast strains are a sister group, i.e. ((b, c),(a, d)). If the JC method is implemented via maximum likelihood (or Bayesian) then it correctly recovers the true tree ((a,b), (c,d)). The artefact is known as long branch attraction.

Clustering algorithms based on a distance matrix that implements the JC correction, tend to perform poorly for long branch-attraction artefacts. This doesn't mean the correction is useless, just not particularly powerful. The issue is that distance matrix methods don't "stick" to the model and the clustering will introduce a layer of inaccuracy. Normally presentation of a "corrected" distance matrix combined with a clustering algorithm will require bootstrapping (resampling with replacement) to assess whether a given cluster is supported. Parameterised distance matrices, using neighbour-joining clustering in in-conjunction with a bootstrap are considered okay.

The @user172818 mentioned parsimony and this method is considered less reliable because it can't implement a JC correction. IMO it is possible weighted parsimony could make a "come-back", but it would be really complicated to implement a biological weighting method and would require extensive, independent calculations.
> 6 votes

# Answer
A great question, though a little ambiguous. I don't know what "general clustering algorithms" refer to. For biological sequences, building a tree can be thought as a way of clustering. Anyway...

There are different tree building algorithms. Max parsimony (MP), max likelihood (ML) and bayesian algorithms directly take sequences as input. They are distinct from distance based clustering.

Then there is a class of distance based algorithms in phylogenetics. They start from an all-pair distance matrix and aim to find a tree that is most compatible with the matrix. Among them, UPGMA is largely hierarchical clustering. Neighbor-joining is somewhat similar to hierarchical clustering, but it builds unrooted trees. FastME uses a very different approach. Briefly, for each topology it tries to find the best branch lengths by weighted least square (for details, see its paper); the best topology minimizes total branch lengths. FastME finds the best topology via nearest neighbor interchange (NNI), closer to some ML algorithms (e.g. PhyML).

> does it brings anything new at the level of clustering algorithm?

Personally, I think the distance based tree building methods are more sophisticated than naive hierarchical clustering, *provided that* the real relationship follows a tree topology.
> 5 votes

---
Tags: sequence-alignment, phylogenetics, phylogeny, clustering
---


# Question
Title: Consensus sequence from SAM or BAM file?
I am trying to perform reference-based assembly. Most of the tutorials teach how to create a bam file and view alignemnts in IGV or Tablet. But, I want a assembled genome sequence in fasta format. How can I do that ? Is there any tool or linux syntax available for this purpose ?

# Answer
> 1 votes
A few years ago I used AlignGraph (paper and code) for a task like this. I gave AlignGraph my raw reads, so I'm not sure if it accepts SAM/BAM input or if you have to go back to the original Fastq reads.

I'm not sure if any related/improved tools have been published since AlignGraph. Internet searches for "reference-guided assembly" return mostly transcriptome assembly tools. It might be worth looking for articles that cite the AlignGraph paper.

---
Tags: sam, genome, assembly
---


# Question
Title: DESeqDataSetFromTximport invalid rownames length
I am trying to use `DESeqDataSetFromTximport` function from `DESeq2` package to construct `dds` object:

```
dds <- DESeqDataSetFromTximport(txi, sampleTable, ~Group)

```

And somehow it is giving me the following error:

> Error in `rownames<-`(`*tmp*`, value = c("ENSMUSG00000000001", "ENSMUSG00000000003", : invalid rownames length

I used the same code for other datasets but for the current one it is somehow not working. My `sampleTable` looks like this:

```
       ID      Group
1    de_1 denervated
2    de_2 denervated
3    de_3 denervated
4    de_4 denervated
5   den_5       sham
6    ho_1        IgG
7    ho_2        IgG
8    ho_3        IgG
9    sh_1       sham
10   sh_2       sham
11   sh_3       sham
12   sh_4       sham
13   sh_5       sham
14 VLA_10  anti-VLA4
15 VLA_12  anti-VLA4
16  VLA_7  anti-VLA4

```

`txi` was formed in the following way:

```
k = keys(txdb, keytype = "GENEID")
df = AnnotationDbi::select(txdb, keys = k, keytype = "GENEID", columns = "TXNAME")
tx2gene = df[, 2:1]  # tx ID, then gene ID
txi <- tximport(names, type = "salmon", tx2gene = tx2gene)

```

Any suggestions would be greatly appreciated

# Answer
I got exactly the same problem not long ago (see my other question what I tried to do). We used the same function!

My code ran fine on my HPC but the same code/data failed on my laptop. I later found out I had to upgrade my R to the latest version (3.5.1) and upgrade all my bioconductor packages to the latest. My HPC was on a recent bioconductor docker, so there was no problem.

Please make sure everything in your environment latest (R, R packages, Bioconductor).
> 3 votes

# Answer
The same error occurred to me with `DESeq2_1.22.2` using `R 3.5.0`. By trial and error, I managed to solve the issue when I updated `SummarizedExperiment` from `1.10.1` to `1.12.0`.
> 2 votes

---
Tags: deseq2
---


# Question
Title: What does "fetching by region is not available for SAM files" mean?
I am used to `gzip`/`biopython` solutions when dealing with sequencing data, but now I wish to switch to more elegant `pysam`. So I looked at the manual, but ran into quite bizarre troubles with the first couple of lines using my bam file

```
import pysam
samfile = pysam.AlignmentFile("3_Tms_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tms_b3v08_scaf000159'):
    print(read)
samfile.close()

```

returns `ValueError: fetching by region is not available for SAM files`. Well, the file is `bam`. I tried to google the error but the only hits I found are the lines in the source code of pysam that check if the file is bam/cram or sam, so somehow `pysam` thinks that my bam is a sam. How can I convince it otherwise? I have also noticed that the manual is for python 2.7, that's maybe where the problem comes from...

# Answer
Your *3\_Tms\_1\_mapped.bam* file, despite its filename extension, is in fact a bgzipped **SAM** file. You can verify this using `htsfile`, which is a small utility packaged with HTSlib:

```
$ htsfile 3_Tms_1_mapped.bam
3_Tms_1_mapped.bam: SAM version 1.3 BGZF-compressed sequence data

```

(For files that really are in BAM format, it reports *BAM version 1 compressed sequence data*.)

So the error message is accurate in this case.
> 5 votes

# Answer
That isn't actually a bam file as John Marshall figured out. I am keeping the rest of my answer since it could be useful to someone else, but the issue here was that you had a compressed (bgzipped) sam file and not an actual bam file and that's why you were getting that error. When I sorted your file in preparation for indexing it, I converted to a bam which is why the rest of my answer worked.

---

You don't have the index file for the bam file you're using. I used this script on the file you linked to (changing the name so that it corresponds to the right file and a contig in that file):

```
#!/usr/bin/env python3

import pysam
samfile = pysam.AlignmentFile("3_Tce_1_mapped.bam", "rb")
for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
    print(read)
samfile.close()

```

The directory I ran it in had:

```
$ ls 3*
3_Tce_1_mapped.bam

```

And I got the error you described:

```
$ foo.py
Traceback (most recent call last):
  File "/home/terdon/scripts/foo.py", line 5, in <module>
    for read in samfile.fetch('3_Tce_b3v08_scaf005149'):
  File "pysam/libcalignmentfile.pyx", line 1107, in pysam.libcalignmentfile.AlignmentFile.fetch
ValueError: fetching by region is not available for SAM files

```

However, indexing the bam file fixed it:

```
$ samtools sort 3_Tce_1_mapped.bam > 3_Tce_1_mapped.sorted.bam
$ mv 3_Tce_1_mapped.sorted.bam 3_Tce_1_mapped.bam
$ samtools index 3_Tce_1_mapped.bam
$ ls 3*
3_Tce_1_mapped.bam  3_Tce_1_mapped.bam.bai

$ foo.py | wc
227    2724   16725

```

So just sort and index your files before attempting to seek in them. Which makes sense since the index's job is primarily to allow seeking.
> 5 votes

---
Tags: python, sam, pysam
---


title: filter secondary alignments using pysam

text:
# Question
Title: filter secondary alignments using pysam
I am new to `python` and trying to learn. The below is an attempt to filter out secondary reads in a `bam` using `pysam`. I ran `mark duplicates` and got an error on several `bam` files due to secondary alignments. I added comments to each line to illustrate my thought process. Thank you :).

The specific `mark duplicates` error:

```
Value was put into PairInfoMap more than once

```

Here is the script:

```
#! /usr/bin/python  ## call python script
import sys ## import python system functions
import pysam  ## import module

bam = pysam.AlignmentFile(*.bam, "rb")  ## open bam and read
click.echo("Reading BAM file")  ## output message
for read in bam:  ## start loop and iterate over each bam
    if read.is_secondary=true  ## not the primary alignment
    if read.has_tag('XS')  ## look for Secondary Alignment Score tag
        XS = read.get_tag('XS')  ## read and store AS value
    if read.has_tag('AS')  ## look for Alignment Score tag
        AS = read.get_tag('AS')  ## read and store value os XS tag
    if AS < XS  ## Alignment score less then
        read.write(read)  ## only print primary alignments
bam.close()  ## end processing

```

# Answer
Quick look at your script suggests that you are comparing the XS and AS from different reads. I would recommend removing the continue from the `read.has_tag('XS')` continue block.

I think that if you filter out `is_secondary` and `is_supplementary` you don't need to do the XS to AS comparison assuming whatever generated for BAM did so correctly.

```
for read in bam.fetch()
    if read.is_secondary=true
        continue
    if read.has_tag('XS') <- This block stores XS and then moves on to the next read
        XS = read.get_tag('XS')
        continue
    if read.has_tag('AS')
        AS = read.get_tag('AS')
    if AS < XS  <- This block compares the AS of the current read to the XS of the previous read that had one.
        read.write(read)
bam.close()

```
> 3 votes

---
Tags: python, pysam
---