# Question
Title: Fix dash/empty ALT alleles for deletions in VCF
I have a VCF file that has deletions specified like this, with the REF allele containing only the deleted portion and the ALT allele just saying `-`:

```
chr1 101 . AAA -

```

I would like to combine this with another VCF in which deletions are specified like this, with REF and ALT both extended by 1bp to the left to ensure neither is empty:

```
chr1 100 . GAAA G

```

Is there a tool to convert from one of these representations to the other? I have already tried `vt normalize` and `bcftools norm` but neither of them seem to do this.

I don't really mind which representation I end up with although I think the second is the one more consistent with the spec.

# Answer
> 6 votes
According to the specs the ALT is not allowed to contain `-`. While technically supported, 'N' as the alternative allele is sometimes used unnecessarily and this script can also fix those cases if one exchanges the two "-" within the fix\_record function for "N".

Here is a little python (\>=3.6) script that is able to fix it. It makes use of the pysam module. You need a reference file and this fasta file must be indexed by `samtools faidx reference.fa`. Please note that it doesn't support multiple values in the ALT column for records that contain a `-`.

**Usage:**

```
python fixvcf.py -f reference.fasta -o out.vcf in.vcf

```

**The script:**

```
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """
        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) -> pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """
        if "-" in record.alts:
            ref = self.fasta.fetch(region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}")
            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see: https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """
        print(self.vcf.header, end="")
        print(*map(self.fix_record, self.vcf), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

```

A lot of the code above is boilerplate for running a script on the command line. The critical code for the transformation is in the `fix_record` function. The relevant code is below, expanded, with line numbers and a description of each command.

```
[1]   ref = self.fasta.fetch(
          region=f"{record.chrom}:{record.pos - 1}-{record.pos - 1}"
      )
[2]   record.pos = record.pos - 1
[3]   record.ref = f"{ref}{record.ref}"
[4]   record.alts = tuple(
          map(
              lambda alt: ref if alt == "-" else f"{ref}{alt}",
              record.alts
          )
      )

```
1. This grabs the base preceding the inserted/deleted sequence.
2. This corrects the position of the variant to account for prepending another nucleotide.
3. This defines the new corrected REF allele: concatenate the base grabbed in \[1\] with the original REF allele
4. Apply an anonymous (lambda) function to each ALT allele; the function will replace `-` with the base grabbed in \[1\], or preprend the base to the allele if it is not `-`

# Answer
> 1 votes
Probably it is better to use dedicated tools, however, here is a quick R workaround. It is relatively straightforward to remove the extra base (first letter of column 4) in the second type of output:

```
# just a workaround to encolose the second file in a data frame
# you can just read your file as a table using fread(), read.csv(), ...
test <- data.frame(
  col1 = "chr1",
  col2 = 100,
  col3 = ".",
  col4 = "GAAA",
  col5 = "G",
  stringsAsFactors = FALSE
)

> test
  col1 col2 col3 col4 col5
1 chr1  100    . GAAA    G

# remove the first letter in column 4
test$col4 <- substring(test$col4, 2)
# converting column 4 to "-"
test$col5 <- "-"

> test
  col1 col2 col3 col4 col5
1 chr1  100    .  AAA    -

```

EDIT:

Here is how it can be done with `awk` for speed and memory concerns: The logic is the same as above, remove the first letter in column 4 and then change column 5 to "-".

```
awk '{$4 = substr($4,2); print}' infile > outfile
awk '{$5 = "-"; print}' outfile > outfile2

```

infile

```
chr1 100 . GAAA G

```

outfile

```
chr1 100 . AAA G

```

outfile2

```
chr1 100 . AAA -

```

# Answer
> 1 votes
One more thing. If empty ALT allele is not "-", you can change it such as "\*", as follows.

```
    if "*" in record.alts:
        ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)
        record.pos = record.pos - 1
        record.ref = f"{ref}{record.ref}"
        record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))

```

# Answer
> 1 votes
This is my first post and I will add some explanations. When I tried to use this script, `print(*map(self.fix_record, self.vcf), sep="", end="")`, it will increase memory usage because it will load everything into memory and then start to apply fix\_record function. So, I just used two lines `for line in self.vcf: print(self.fix_record(line), sep="", end="")`. Because `pysam` uses generator by default, if we process one by one by pysam's `VariantFile` function, it will consume very little memory. In my case, the VCF file was more than 200GB, so I had to fix this script for my VCF file.

Second my VCF file used "\*" for empty ALT, so I fixed two lines (1) `if "*" in record.alts:` (2) `record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))`

```
#!/usr/bin/env python
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """

        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) ->     pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """

        if "*" in record.alts:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "*" else f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see:
https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """

        print(self.vcf.header, end="")
        #print(*map(self.fix_record, self.vcf), sep="", end="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

```

# Answer
> 0 votes
Edited version of Namshin's code this can handle if REF or ALT is '-' ( these records in vcf file are observed when there is insertion/deletion). Added documentation in the script to explain code block I have modified. Hope this useful. Tested example with GRCh37

```
python fix_vcf.py -f genome.fa -o out_fixed.vcf input.vcf

```

example input.vcf:

```
##fileformat=VCFv4.2
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . - T 1 PASS  NS=1  GQ  10
9 98218616  . - C 1 PASS  NS=1  GQ  10
17  41245200  . AT  - 1 PASS  NS=1  GQ  10
13  32929213  . C - 1 PASS  NS=1  GQ  10

```

out\_fixed.vcf:

```
##fileformat=VCFv4.2
##FILTER=<ID=PASS,Description="All filters passed">
#CHROM  POS ID  REF ALT QUAL  FILTER  INFO
7 6027157 . T TT  1 PASS  NS=1
9 98218616  . A AC  1 PASS  NS=1
17  41245199  . GAT G 1 PASS  NS=1
13  32929212  . CC  C 1 PASS  NS=1

```

```
# -*- coding: utf-8 -*-
import argparse
import sys

import pysam

class Application:
    def __init__(self):
        self.args = Application.get_args()
        self.vcf = pysam.VariantFile(sys.stdin)
        self.fasta = pysam.FastaFile(self.args.fasta)

    @staticmethod
    def get_args() -> argparse.Namespace:
        """
        Read command line arguments and set input and output sources.

        :return: parsed command line arguments
        """

        parser = argparse.ArgumentParser(prog="fixvcf.py", description="Fixes deletion represented by -")
        parser.add_argument('--version', action='version', version='%(prog)s 0.1')

        parser.add_argument("input", help="input file, use - to read from stdin")
        parser.add_argument("-f", "--fasta", help="fasta reference file", required=True)
        parser.add_argument("-o", "--output", help="output file")

        if len(sys.argv) == 1:
            parser.print_help(sys.stderr)
            sys.exit(1)

        args = parser.parse_args()

        if args.output:
            sys.stdout = open(args.output, "w")

        if args.input != "-":
            sys.stdin = open(args.input, "r")

        return args

    def fix_record(self, record: pysam.libcbcf.VariantRecord) ->     pysam.libcbcf.VariantRecord:
        """
        If the ALT column of the records contains a `-`, the variant position is shift one position to the left.
        For reaching this the reference base for the position is determined and prepended to the current REF bases.
        The ALT value `-` is replaced with the determined ref base on position before. In case there are multiple ALT
        values the determined ref base is prepended to them.

        If the REF column of the records contains a `-`, the variant position will change.
        For reaching this the reference base for the position is determined and prepended to the current ALT bases.
        The REF value `-` is replaced with the determined ref base on position. In case there are multiple ALT
        values the determined ref base is prepended to them.

        :param record: vcf file record
        :return: fixed vcf record
        """

        if "-" in record.alts:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 2, end=record.pos -1)

            record.pos = record.pos - 1
            record.ref = f"{ref}{record.ref}"
            record.alts = tuple(map(lambda alt: ref if alt == "-" else f"{ref}{alt}", record.alts))

        if "-" in record.ref:
            ref = self.fasta.fetch(reference=record.chrom, start=record.pos - 1, end=record.pos)

            record.pos = record.pos
            record.ref = f"{ref}"
            record.alts = tuple(map(lambda alt: f"{ref}{alt}", record.alts))

        return record

    def start(self):
        """
        Print the header of the vcf file. Afterwards the fix_record method is applied to each record in the vcf file.

        For more information about python's map() see:
https://docs.python.org/3.6/library/functions.html#map

        For more information about the unpacking mechanism using `*`see:
        https://docs.python.org/3/tutorial/controlflow.html#tut-unpacking-arguments
        """

        print(self.vcf.header, end="")
        #print(*map(self.fix_record, self.vcf), sep="", end="")
        for line in self.vcf:
            print(self.fix_record(line), sep="", end="")

if __name__ == "__main__":
    app = Application()
    app.start()

```

---
Tags: vcf
---


# Question
Title: Fastq: how can I check if they are from DNA or RNAseq data?
I have (gave me) Illumina fastq files, which I want to use for variant calling, and I do not know if they are DNAseq or RNAseq data. How can I check this? I do not have any report or who to ask.

Many thanks

# Answer
There are a few quick'n'dirty ways depending on the type of data. In any case you want to align your files to a reference genome and then check the distribution of reads, either on a genome browser or with tools such as RSEQC which calcualtes the fraction of reads aligning to exon, intron, intergenic etc.

RNA-seq, if you use a standard aligner such as bowtie2 or bwa, and you have a higher species that splices its RNA, then you should see most reads aligning to exons, and quite a fraction of reads that span a splice junction being unmapped (because the exon-intron-exon gap is large and non splice-aware tools do not "bridge" this without tremendous drop in mapping quality so the reads goes unmapped).

DNA-seq (if whole genome) should have a somewhat even coverage across the genome on the global scale. Targeted assays such as ATAC-seq should give distinct peaks that you can easily see in a genome browser just by eye.

If it is exome-seq (so targeted DNA-seq via exon pulldown) then it might get tricky because also most reads (like in RNA-seq) should align to exons, but you should not have spliced reads so the 5' part of the read starting in exon1 and 3' part ending in exon2. Instead many reads should also cover the exon-intron boundary.

As said, I would simply align it and check on a browser, I could imagine one can make a sophisticated guess by just looking at the alignments by eye.

A different topic is reliability, because I personally would not touch a dataset with unknown origin, because that also means you have no idea about the protocols that were used for library prep and what the sample actually is, so I wonder which question your analysis is going to answer.
> 7 votes

# Answer
Honestly, ask who made them. If the reads are from an organism without splicing, it will not be trivial to figure it out with 100% certainty.

If you are calling variants, I don't think it matters much, except that higher organisms have such a small % of their genome transcribed, it should be obvious very quickly if you have the desired coverage across the entire genome, or just the transcripts.
> 4 votes

---
Tags: rna-seq, phylogenetics, fastq
---


title: Build tree to calculate beta diversity

text:
# Question
Title: Build tree to calculate beta diversity
I am trying to calculate beta diversity\* using the UniFrac distance on my data, each sample having OTUs values.

To calculate it, I saw it required building a tree with percentages that represented the OTUs. However, I was unable to figure what percentage fit each OTU to build the tree from the instructions I found on the internet.

For example, I saw in this link: http://readiab.org/book/0.1.1/3/1#4.1.2, that they build this tree:

Could you please explain to me how to build the tree?

---

* This is the diversity between local sites and enables the total in a landscape to be calculated (gamma). The tree would therefore appear (presumably) to represent the diversity between sites.

# Answer
> 0 votes
Can you give a little more information about what tools and data sets you are using? How did you generate the OTU table etc?

If the underlying data is nucleotide data and this is a phylogentic tree. It may well be a population genetic tree based for example in allele frequencies.

For classical phylogenetics then FastTree is a very simple tool to generate a decent phylogeny. If you have the sequences of each OTU, you can align them in MegaX with muscle, or simply using muscle from the command line if you prefer. You can then generate the tree with this command (presuming you are on windows):

```
FastTree.exe -gtr -nt < aligned_sequences.fasta > otus_tree.tre

```

I hope that helps!

---
Tags: phylogeny, phylogenetics
---


# Question
Title: Will using smaller kmers help get larger contigs? If not, then what?
I've been using MEGAHIT to assemble metagenomes, with particular focus on specific genomic areas.

Sometimes all I get is gene fragments or pathway fragments (eg. if I know that genes A, B, C, D and E should be together, I only get A, B, C in one contig, and maybe D and E in another). That is understandable when the sequencing depth is low, and I'm assume that it's the best I can get.

However. In two of my metagenomes, the genes of interest seem to be abundant (depth of 100-300+ listed for MEGAHIT contigs). And I still have the same problem. I assume that these 'breaks' between contigs are due to natural sequence variation which breaks the assembly process. (while sequencing errors are also possible, I've used bbduk to clean my data prior to assembly)

I'd like to have larger contigs because I'm curious about the gene order. Things I've tried (did not work): visualizing assembly graphs with Bandage around Blast hits, starting with a smaller k-min (21 as opposed to 27).

So I guess I'm wondering if using even smaller k-mers might make a difference. Or if someone has another suggestion, I'm happy to listen! Thank you.

Megahit setting used:

```
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 300 --presets meta-large -o output
megahit -r input_files.fastq --num-cpu-threads 32 --min-contig-len 1500 --presets meta-sensitive -o output

```

Meta-large starts from minimum kmer size 27; meta-sensitive from 21.

# Answer
You can try s-aligner (free for 15 days) which usually gets quite larger contigs and quite larger NG50 for metagenomic data containing viruses (also phages). Indeed metaSPAdes also gets larger contigs, despite shorter than s-aligner.

Ref. s-aligner: a greedy algorithm for non-greedy de novo genome assembly
> 0 votes

---
Tags: gene, assembly, pathway, k-mer
---


# Question
Title: Explanation for RNA-seq samples not clustering in PCA as expected
A colleague is analysing RNA-seq data - the study design is 2 treatments, 3 replicates, 3 tissues. In their PCA plot the samples clustered neatly by tissue. Except for two samples - two tissue samples originating from the same animal clustered with the wrong tissue.

I suggested that the two samples had likely been switched accidently somewhere along the way. They think its explained simply by biological variation which I'm not convinced of. Are there any technical explanations for why samples could cluster on a PCA in such a way?

# Answer
> 1 votes
It's hard to say what is going on without looking at the plots. But it is a reasonable guess that samples are mislabelled and I think it is worth checking the experimental log carefully. It is also possible that the two samples are of low quality. People often use PCA to check and remove outliers from their RNA-seq libraries. So you may also check the RNA-seq library quality for those two samples.

# Answer
> 1 votes
1. **Check for mislabelled samples**

For the mismatch you can check if there is a marker gene for these 3 samples for that particular cohort/tissue that you actually think should be expressed in this tissue only and not in other cohort.

2. **Check for outliers**

To check the logs of mapping reads of the samples that aren't clustering within expected tissue it might be due to low reads. You can use hierarchical clustering and project it in the form of a heatmap. Check if these two samples have varying expression levels to its corresponding biological replicates.

# Answer
> 1 votes
Have you corrected for batch effects? This is usually a reason, as it can mean that the variation being observed from your PCA is due to other factors rather than tissue and/or treatment as you expect - of course, another can be they are incorrectly labeled !

If the samples were sequenced at different times, or by different services, this can often happen. If that is the case, what you will usually find are samples being clustered based on the sequencing service, for example.

Use `ComBat` or `limma` packages in R to fix variations in the data that may be due to this, and rerun the PCA. Best of luck.

---
Tags: rna-seq, pca
---


# Question
Title: Question about public availability of human SNP dataset with country of origin
Is there a publicly available data set for humans SNPs, preferably together with the country of origin? If so, could someone please point me towards it.

I've looked through other questions on the site and have been unable to find a similar one. I don't know how sensitive human SNP data is, so I'm not sure if this would be public.

# Answer
> 3 votes
1000 genomes, HGDP and SGDP all contain SNP data from individuals where the country of origin is known.

They all have slightly different numbers of people from different populations (1000 genomes and HGDP have more individuals in each pop, but fewer pops, and SGDP is the opposite).

---
Tags: database, snp, public-databases, human-genome
---


# Question
Title: Converting FASTA Sequences to SNP FASTA Format
I have two FASTA sequences. Each is from a different species' mitochondrial genome. The sequences are pasted below.

```
>speciesA_control_region_partial_sequence_mitochondrial
TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA
>speciesB_control_region_partial_sequence_mitochondrial
TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA

```

My goal is to extract all SNPs in these sequences and format them in SNP FASTA format, as shown below:

```
>rs11540652
GAACC[A/G]GAGGCCC

```

I specifically need this format (rather than, say, VCF) so I can search for RFLP sites using the following website: http://bio.kuas.edu.tw/snp-rflping2/rflpUI.jsp

Does anyone know of an existing program which will accomplish this? Also, what additional materials do I need to detect these SNPs? I presume that I need assembled mitochondrial reference genomes for each species - anything else?

# Answer
Your question is not very clear. Where does the flanking sequence in :

```
>rs11540652
GAACC[A/G]GAGGCCC

```

come from? I'll assume its meant to be the reference sequence and that speciesA is the reference and that the first nuc in the \[\] is from the reference... I don't know an existing existing program but here's a bespoke one.

```
speciesA='TTAAAGTGATGTAATCCTACATTAATGAAAAATCAAAAGTTATAAGAACTTAAATACATTATATCATCAAATAAATATGAAGGTAGACATAAACCACTAAACTTCAAACTCCATTAAATATGTTATAAAAATGACGATATTGAATTGTCCTATCATAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTCAGAGTCAAATGTAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCAAGGACAAAAATCGTGGGGGTAGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCCTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCATTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATTAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATTATTTTTTCCCCTGACTTTCCTGTCAAGCACCCCTCTCGGCTTTTACGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
speciesB='TTAAAGTAATGTAACCCTACATTAATGAAAAATCAAAA-TTATAGGAACTTAAATACATTACACCATCAAATAAATATGAAGGTAGACATAAACCATTGAATTTTAAATTTCATTAAACATGTTAAAAAAATGACGATATTGAATTGCCCTATCACAACTCTCATCAGTCTAGATATACCAGGACTCACCACCTCTGCAAGTAAGAGTCAAATGCAGTAAGAGACCACCATCAGTTGATTTCTTAATGTACACGTTTATTGATGGTCACGGACAAAAATCGTGGGGGTCGCACACAGTGAACTATTCCTGGCATTTGGTTCCTATTTCAGGTTCATCAATTGAATACATTCCCCTAACATTTCTTGACGCTTGCATAAGTTAATGGTGGAGTACATACTCCTCGTTACCCACCATGCCGAGCCTTCACTCTATAGGGCTACTGGTATCTTTTTTCTATTTCCTTTCACTTGACATTTCAGAGTGCATACAGAAATGTTATATCAAGGTTGAACATTTCCTTGCCCGGCGGAAAATGTATGAATGTAATTAAACTTTTATTTATGAATTGCATAACTGATATCAAGAGCATAAAGTATCATTTTTTCCCCTGACTTTCCTATCAAGCACCCCTCTCGGCTTTTGCGGGCCAAAACCCCCCCTCCCCCCTAAACTCGTGAGATTCCATTGTGTCTGCAAACCCCCCGGAAACAGAGCAAATCTTACTAGTTTTATCCCTTCTCAAATTTTTGTTCACTTGTATTTTTA'
counter=0
for i, nuc in enumerate(speciesA):
    if nuc != speciesB[i]:
        counter+=1
        print ('>rs'+str(counter))
        print (f'{speciesA[i-5:i]}[{nuc}/{speciesB[i]}]{speciesA[i+1:i+6]}')

```

Gives:

```
>rs1
AAAGT[G/A]ATGTA
>rs2
TGTAA[T/C]CCTAC
>rs3
CAAAA[G/-]TTATA
>rs4
TTATA[A/G]GAACT
>rs5
CATTA[T/C]ATCAT
>rs6
TTATA[T/C]CATCA
>rs7
AACCA[C/T]TAAAC
>rs8
CCACT[A/G]AACTT
>rs9
CTAAA[C/T]TTCAA
>rs10
AACTT[C/T]AAACT
>rs11
TCAAA[C/T]TCCAT
...

```

As for 'I presume that I need assembled mitochondrial reference genomes for each species - anything else?' These sequences have clearly been aligned. So there's that.

I hope that helps.
> 2 votes

---
Tags: fasta, snp, genomics
---


# Question
Title: Confirmation about PERMANOVA and Variance partitioning (permutest)
Could you please help me to understand what is the difference between Permutational multivariate analysis of variance (PERMANOVA,adonis) and Variance partitioning (permutest).

I understood that permutest is used to confirm the PERMANOVA result that differences in PERMANOVA is not due to dispersion.

```
###adonis_calculation
BC <- phyloseq::distance(phyloseq_prop, "bray")
pmv =adonis(BC ~Season*Region, as(sample_data(physeqN3), "data.frame"),permutations = 5000)
pmv

Call:
adonis(formula = BC ~ Season * Region, data = as(sample_data(physeqN3),      "data.frame"), permutations = 5000)

Permutation: free
Number of permutations: 5000

Terms added sequenti
Total         99   19.8578                 1.00000 ally (first to last)

              Df SumsOfSqs MeanSqs F.Model      R2 Pr(>F)
Season         1    1.8586 1.85862  12.036 0.09360  2e-04 ***
Region         1    2.4218 2.42183  15.682 0.12196  2e-04 ***
Season:Region  1    0.7522 0.75222   4.871 0.03788  2e-04 ***
Residuals     96   14.8251 0.15443         0.74656

```

Many thanks

# Answer
Permutational Multivariate Analysis of Variance (PERMANOVA) assumes no distribution which means it is a non-parameteric statistical test. According to the assumption of PERMANOVA test, each group must have same dispersion in the data. Before performing any statistical test, we must check the underlying assumptions. Variance partitioning (permutest) is a confirmatory test to check the assumption (similar dispersion). You can use permutest.betadisper for permutest. After that you can perform PERMANOVA test using adonis().

Interpretation of PERMANOVA test –

1. If p value is small (\<0.05 ) we reject the null hypothesis (No difference between the groups) otherwise we accept the null hypothesis (Significant difference between the groups).
2. Coefficient of determination, R-squared value indicates the variation in the distances being explained by groups under the study. Higher value indicates a better fit.

You can further look into the confidence interval for each group.

Hope this confirms the answer to your question.
> 1 votes

---
Tags: r, statistics, phyloseq
---


# Question
Title: How to grep with a list of genes from a gff the GO terms
I have a text file with a list of name of genes (alias names) like this (`**` added for emphasis):

```
Solyc00g019980
**Solyc01g005020**
..
.

```

And i have a gff file like this:

```
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=**Solyc01g005020**;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;**Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085**
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

```

I would like to grab only the ones i have in my list gene text file (Alias names) and the corresponding "Ontology\_term" next to it.

Terminal is more convenient if there is a solution like grep or something.

Thanks in advance

edit: i would like to get this output, as an example

```
Solyc01g005070  Ontology_term=GO:0005543
Solyc01g005230  Ontology_term=GO:0008168

```

i dont mind if this is included "Ontology\_term=", i would prefer not to be included though

thanks again

`grep -f list.txt data.gff` doesn't work for me, because it grabs all the other as well.

# Answer
@Jim Hat

I took a stab at this given the example files you included.

If the 'list.txt' file looks like this:

```
Solyc00g019980
Solyc01g005020

```

And the data.gff file looks like this:

```
SL4.0ch01   maker_ITAG  gene    46201   59605   .   -   .   ID=gene:Solyc01g005020.3;Alias=Solyc01g005020;Name=Solyc01g005020.3;length=13404
SL4.0ch01   maker_ITAG  mRNA    46201   59605   .   -   .   ID=mRNA:Solyc01g005020.3.1;Parent=gene:Solyc01g005020.3;Name=Solyc01g005020.3.1;Note=Salt Overly Sensitive 1;_AED=0.25;_QI=241|1|1|1|0|0|23|311|1151;_eAED=0.25;Ontology_term=GO:0006812,GO:0015299,GO:0016021,GO:0055085
SL4.0ch01   maker_ITAG  three_prime_UTR 46201   46511   .   -   .   ID=three_prime_UTR:Solyc01g005020.3.1.0;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    46201   46844   .   -   .   ID=exon:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 46512   46844   .   -   0   ID=CDS:Solyc01g005020.3.1.1;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    47281   47585   .   -   .   ID=exon:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  CDS 47281   47585   .   -   2   ID=CDS:Solyc01g005020.3.1.2;Parent=mRNA:Solyc01g005020.3.1
SL4.0ch01   maker_ITAG  exon    48276   48543   .   -   .   ID=exon:Solyc01g005020.3.1.3;Parent=mRNA:Solyc01g005020.3.1

```

(Note I removed the \**s, which I assumed were added for emphasis)

You can do the following:

```
aliasNames=$(grep -f list.txt data.gff | cut -d':' -f2 | cut -d';' -f1)
ontologyTerms=$(grep -f list.txt data.gff | cut -d';' -f 8 | cut -d'=' -f2)
paste <(echo "$aliasNames") <(echo "$ontologyTerms")

```

The first line filters the file based on your list and then extracts the alias names from the filtered GFF file. The second line filters the GFF and extracts the ontology terms (after the equal sign). And the third line joins the two as columns.

The result that is printed to the screen is:

```
Solyc01g005020.3
Solyc01g005020.3.1      GO:0006812,GO:0015299,GO:0016021,GO:0055085
Solyc01g005020.3.1.0
Solyc01g005020.3.1.1
Solyc01g005020.3.1.1
Solyc01g005020.3.1.2
Solyc01g005020.3.1.2
Solyc01g005020.3.1.3

```

Given the examples you included all of the lines of the GFF are included, but not all have gene ontology information.
> 2 votes

---
Tags: grep
---


# Question
Title: Extract features from fasta sequences and train the classifier
I am new to the bioinformatics field. I have positive and negative protein sequences for acetylation PTM. Now, I want to train a classifier, say SVM. What will be the next step? How can I convert these sequences into usable features? Any information or links would help.

```
>P31327_55|1|testing
KAQTAHIVLEDGTKMKGYSFGHPSSVA
>P31327_57|1|testing
QTAHIVLEDGTKMKGYSFGHPSSVAGE
>P31327_119|1|testing
APDTTALDELGLSKYLESNGIKVSGLL
>P31327_157|1|testing
LATKSLGQWLQEEKVPAIYGVDTRMLT

```

# Answer
> 1 votes
Just a recommendation: Be careful when making your data sets.

P31327\_57, P31327\_119, P31327\_157 are all annotated as "By similarity" to the mouse entry in Swiss-Prot https://www.uniprot.org/uniprot/P31327#ptm\_processing, which means there is no direct experimental evidence in the human entry. If you do have publications that report experimental evidence for these PTMs, please don't hesitate to share them with the UniProt team.

I also recommend reading this UniProt help page about using UniProt data in negative datasets if you haven't already done so: https://www.uniprot.org/help/negative\_datasets

---
Tags: fasta, proteins, machine-learning, features
---


# Question
Title: In SBOL, is there a way to distinguish biological replicates?
I have a library of biological implementations (i.e. plasmids) of genetic circuits and their experimental data associated. I want to implement an *in silico* version with SBOL, but to be consistent with what I have I'm my freezer, I need to identify not only the genetic circuit but also the biological replicate (i.e. the different colonies picked), as they have their own associated experimental data.

# Answer
> 6 votes
Support for distinguishing replicates was added in SBOL 2.2, along with other support for linking data throughout the experiment cycle. I will explain this in terms of SBOL 3, which is simpler, then give the SBOL 2 translation at the end.

For SBOL 3, all designs, including those for a plasmid, a strain, etc, are specified in terms of the `Component` class. This is just a "paper" design, however, and may or may not actually exist. When you have a real physical aliquot that you want to talk about, this is represented by an `Implementation` object: this is your biological replicate.

The `Implementation` doesn't really have any content of its own; its nature is defined by linking back to the intended design with a `prov:wasDerivedFrom` link. If you've verified that what you got is what you intended, you can supplement that with a `built` link as well.

One additional note, however: when you're representing a colony, you're probably not actually representing just a plasmid, but rather are representing a particular strain that has been transformed by a plasmid. The actual data you have may then be one level further removed, since the data has been gathered from a sample under particular growth conditions. SBOL can represent that as well, by describing the sample as a mixture of particular media, inducers, etc., with a transformed strain, which is itself a cell type containing your plasmid.

Here is an example diagram showing just this sort of relationship, of data associated with three replicates of a sample.

You can find this example in the SBOL3 tutorial from COMBINE 2020.

Finally, in SBOL 2, this works exactly the same way, except that instead of a `Component` your implementation points to a `ComponentDefinition` or a `ModuleDefinition`, depending on the nature of the design you're representing. I won't go into that distinction, because it was a bad idea and there's a reason we got rid of it for SBOL 3.

---
Tags: sbol
---


title: Doubt about using TPM for statistics

text:
# Question
Title: Doubt about using TPM for statistics
We designed an experiment to explore the potential role of carbon dioxide on algae physiology using RNA-Seq. We analyse the differential gene expression using DESeq2 but now we are interested into analyse a few genes in detail.

My PI say to make the statistics over the ration response (expression\_treatment/expression\_control - 1) and my question is the following: Can I use TPM to calculate this ratio and make following statistics? I read in the documentation of DESeq2 that TPM we can't use this normalisation method to differential expression and my my other question is Why shouldn't?

My advisor ask to re-do the statistical analysis using the ratio response (rr) so i need the rr for each sample. I think this doesnt make sense cause the result will be or should be the same that using DESeq2.

Update: I'm using the fold-change calculated by DESeq2 and I'm not sure if I have transformed the data correctly. I took the two columns `log2FoldChange` and `lfcSE` (both are log-transformed variables) so I went back to the non-transformed variable as follow: `FC=2^log2FoldChange` and `SE=lfcSE*2^log2FoldChange`. I saw this approximation in this post. Is that correct? I have obtained a very broad range of error/

# Answer
This is a common question, and the answer is while you *can* calculate this, it's not statistically robust. You're likely to arrive at false conclusions.

Why is this? Because while TPM is something you can calculate from a single sample, your question involves multiple replicates and those replicates vary between each other. Counts in units of TPM are missing information about the variation between samples. This extra information affects fold change estimates, and is why differential analysis of sequencing data is complicated.

See this answer to a related question, and this often-cited blog post for more information.
> 4 votes

---
Tags: rna-seq, normalization, tpm, fold-change
---


title: add task to nextflow and resume

text:
# Question
Title: add task to nextflow and resume
in test.nf (dsl2) I have

```
   ...
   workflow{
       a = Channel...
       processA(a)
       processB(processA.out.emitA)
       processC(processB.out.emitB)
    }

```

The pipeline completed successfully.

Now I want to add `processC_2` that takes the output of processB.

How do I run the workflow to bypass processA execution?

I added processC\_2 to test.nf and then did:

```
nextflow run test.nf -resume run-name

```

But it doesn't run the newly add process.

How do I achieve my goal?

Thank you!

# Answer
It is actually very hard to answer this using only the information you have provided. It'd be necessary to actually see your `test.nf` file and read the processes. The execution seems right and I'd guess that the processB may be not emiting any output as shown in https://www.nextflow.io/docs/latest/process.html#outputs.

It is also possible that when using the `-resume` together with a specific run-name you may be recovering the state of that specific run as shown here https://www.nextflow.io/blog/2019/demystifying-nextflow-resume.html. You can try to execute it without the run-name or without the `-resume` flag so you can check whether the problem is with the emition and channelling of outputs between processes, or with the `-resume` flag.
> 3 votes

---
Tags: nextflow
---


# Question
Title: How to interpret track height in Chipseq?
I recently received some .bigwig files from a chipseq experiment. I loaded the files into IGV.

The interpretation is rather intuitive.

My question is the following: What exactly do the track heights represent?

For IgG control, the track height is low ranging from 0~5.

For positive controls, such as H3K4me3, I see clear peaks with height ~100.

Do track heights simply mean the number of reads mapped to that genomic interval during peak calling? Does it literally mean there are ~100 reads that aligned to that region in the H3K4me3 .bam file?

Thanks so much.

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

Does this make sense to you?

---
Tags: chip-seq, peak-calling
---


# Question
Title: How do you find differences when your metagenome sample size is small?
I received a matrix of genomics features (KEGG annotated metagenomes) abundance from 6 samples belonging to 2 groups. My collaborator is interested in finding differences between the two groups.

The matrix, after scaling and removing zero-variance, is roughly 4000 features. Via shapiro test I already know some 30% of the features are not normally distributed across the samples, so I assume PCA is out of option and any test with the normal distribution assumption.

I already made some preliminary analysis using NMDS to cluster my samples and was thinking of trying also t-SNE and recursive feature elimination. However I am not sure how to find "significant" differences given the small sample size. Any suggestion? Would something like LEfSe work in this instance?

# Answer
> 0 votes
I think there are several ways to figure this out. I wouldn't worry too much about normality as most of the methods are either nonparametric or pretty robust to violations.

1. Since you are using KEGG, I would suggest working at the pathway or module level of annotation (if your features are e.g. KEGG Orthology groups). In principle, you can then combine information across features as well as across samples. Simplest I can think of, without even looking at the literature, would be to use permutations (shuffle row and column labels) to look for modules different in abundance between groups. (Many existing methods are doing some more sophisticated version of this.)
2. I think that LEfSe in theory should work quite well here, it's been a while since I thought about it but I believe this is the problem that it is designed to solve, assuming that I understand your question.
3. There are a lot of methods that are previously described for doing this kind of work. I would suggest reading some reviews of methods, for example the section "Multivariate differential abundance testing" in this article, or here, looking for some inspiration. To some extent it depends on which questions you are interested in answering and what exactly your features represent.

---
Tags: statistics, genomics
---


# Question
Title: Biopython BioSQL error : AttributeError: 'DBSeq' object has no attribute '_data'
Hi I am playing with Biopython and BioSQL,

I manage to create an MySQL database using Managing local biological databases with the BioSQL module tutorial:

and load 2909 fasta sequences in it through:

```
db = server["test_1"]

count = db.load(SeqIO.parse("Newtextfile-03.txt", "fasta"))
        print("Loaded %i records" % count)
        server.commit()

```

from this file: Newtextfile-03.txt

problem is trying just to print the database index list with:

```
for i in db:
    print(i)

```

works , while trying print element by element:

```
for i in db:
    print(i, db[i])

```

or with:

```
for key, record in db.items():
            print ("Key %r maps to a sequence record  %s" % (key, record))

```

gives an error:

```
2557 Traceback (most recent call last):
  File "./try_001.py", line 142, in <module>
    f()
  File "./try_001.py", line 120, in f
    print(i, db[i])
  File "......../lib/python3.8/site-packages/Bio/SeqRecord.py", line 652, in __str__
    lines.append(repr(self.seq))
  File "........//lib/python3.8/site-packages/Bio/Seq.py", line 109, in __repr__
    return f"{self.__class__.__name__}({self._data!r})"

AttributeError: 'DBSeq' object has no attribute '_data'

```

the error occurs for three records namely : \[2164, 2522, 2557\]

inspecting the Newtextfile-03.txt file results in finding fasta sequences of less than 100 aminoacids.

Is there any way to resolve the error ? It sounds strange to me that I can load the records into the database but not retrieving them. Is anything I am missing about the record load ? or similar ?

I am working with the database using:

```
 server = BioSeqDatabase.open_database(
    driver='MySQLdb' ,
    user=password.userz,
    password=password.pazzword,
    host="XXXXXXXX",
    db="bio-data-first")

db = server.new_database("test_1", description="Just for testing")
        server.commit()

```

using driver="MySQLdb", wasn't able to have driver='mysql.connector' to work,

biopython is version 1.78

# Answer
**It's a bug in Biopython 1.78.**

What those three sequences have in common is that they're less than 61 characters long. When you print the record it implicitly calls `__str__` for the record, which calls `__repr__` for the Seq. DBSeq gets its method from Seq, which looks like:

```
    def __repr__(self):
        """Return (truncated) representation of the sequence for debugging."""
        if len(self) > 60:
            # Shows the last three letters as it is often useful to see if
            # there is a stop codon at the end of a sequence.
            # Note total length is 54+3+3=60
            return f"{self.__class__.__name__}('{str(self[:54])}...{str(self[-3:])}')"
        else:
            return f"{self.__class__.__name__}({self._data!r})"

```

...and then it only tries to access the nonexistent `_data` for the short sequences.

So there's a mismatch between DBSeq's lack of a `_data` attribute and its use of `Seq.__repr__`. I don't see a DBSeq class at all in the latest Biopython, so maybe this is fixed already. But there's no tag for 1.79 yet that I can see and I can't get the very latest code to import successfully so I'm not sure. In the meantime if you avoid letting `repr(record)` happen you should be able to work around it.
> 1 votes

---
Tags: fasta, python, database
---


# Question
Title: Common gene naming conventions and how to convert between them
I have two questions, one is a direct question about some RNA expression data. The other is more broad, but motivated by the first. I have downloaded an expression set for a class project and I am unfamiliar with the naming convention, I recognize entrez ID but not this one.

```
    > rownames(Rao)[1:10]
 [1] "0610007P14Rik" "0610009B22Rik" "0610009D07Rik" "0610009L18Rik" "0610009O20Rik" "0610010B08Rik" "0610010F05Rik"
 [8] "0610010K14Rik" "0610011F06Rik" "0610012G03Rik"

```

So my first question is what is this gene naming convention. The second is what are many other conventions I should recognize and what tools can I use for converting between them? Lastly what kind of errors may crop up when I do that and how can I check for mismappings or non-unique mappings?

I am primarily interested in R but methods in python would also be appreciated! Its always good to have more tools.

# Answer
> 1 votes
Biologists are familiar with `SYMBOL` (for good reasons).

Bioinformaticians like `ENTREZID` `ENSEMBL` `REFSEQ` `GENBANK` etc, I'm definitely missing some here.

You don't need to know these by heart, just know that databases don't map 1:1.

I've worked with mouse RNAseq data before, so I agree with user::RamRS comment that it looks like a mouse gene. What type of identifier is it? I had no idea.

You can find this out at NCBI::GENE, by querying for, let's say `Rao[1]` `"0610007P14Rik"`. (https://www.ncbi.nlm.nih.gov/gene/?term=0610007P14Rik)

NCBI::GENE tells you the `SYMBOL::Erg28`, `ENTREZID::8520`, and in the far right a column called `ALIAS::0610007P14Rik`.

Now that you know your identifier type (ALIAS), you can map between databases in R. There are many packages to do this. Someone mentioned biomart in the comment.

I use `AnnotationDbi`

```
library(AnnotationDbi)
library(org.Mm.eg.db)

Rao <- c("0610007P14Rik", "0610009B22Rik", "0610009D07Rik", "0610009L18Rik", "0610009O20Rik", "0610010B08Rik", "0610010F05Rik", "0610010K14Rik", "0610011F06Rik", "0610012G03Rik")

#use the function select to retrieve the column type "SYMBOL" with keys "RAO".
#"RAO" is a keytype of "ALIAS". And I want to search in the mice database "org.Mm.eg.db"

AnnotationDbi::select(x = org.Mm.eg.db,
+                       keys = Rao,
+                       keytype = "ALIAS",
+                       columns = "SYMBOL")
'select()' returned 1:1 mapping between keys and columns
           ALIAS        SYMBOL
1  0610007P14Rik         Erg28
2  0610009B22Rik 0610009B22Rik
3  0610009D07Rik         Sf3b6
4  0610009L18Rik 0610009L18Rik
5  0610009O20Rik         Dele1
6  0610010B08Rik        Gm2026
7  0610010F05Rik 0610010F05Rik
8  0610010K14Rik 0610010K14Rik
9  0610011F06Rik       Mettl26
10 0610012G03Rik 0610012G03Rik

```

In the cases when they don't map 1:1, the function will throw a warning, and you can access it by `warnings()`. It is rarely important though.

Hope this answers the question.

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


# Question
Title: How to check what software has been used to generate bam file?
I have several bam files from unknown origin. I'd like to know what software was used to obtain those bam files, e.g. `bwa-mem`. Is it possible to somehow check it?

# Answer
If you're lucky the information is stored in the BAM header in the lines starting with `@PG` which stores the command lines run on the BAM, use `samtools view -H your.bam` to display the header.
> 6 votes

---
Tags: sam, ngs, read-mapping
---


# Question
Title: Nextflow error: can not invoke method minus () on null object
Hi I was trying to use plink in nextflow and getting the error mentioned in the subject. Here is my code, any help what I am missing. I am new in nextflow, it was ok when input was a single file, but for grouped files, I am getting lost.

```
params.dir='data/'
params.publish="output/"
Channel
.fromFilePairs("${params.dir}*.{bed,fam,bim}",size:3){file -> file.baseName}
.set {plink-data}

process test1{
publishDir params.publish
input:
set t1,file(plk1) from plink-data
output:
file "${t1}.{bed,fam,bim}" into test1_results
script:
output_file = "${t1}"
base = plk1[0].baseName
"""
module load plink
plink --bfile $t1 --maf 0.05 --make-bed --out ${output_file}
"""
}
process test2{
publishDir params.publish
input:
set file(bed), file(fam), file(bim) from test1_results.collect()
output:
set t2,file '${plk2}.{bed,fam,bim}' into test2_results
script:
output_file = "{t2}"
base = plk2[0].baseName
"""
plink --bfile $t2 --maf 0.01 --make-bed --out ${output_file}
"""
}

```

# Answer
Channel names must consist only alphanumeric characters including underscores. Your channel called 'plink-data' contains a hyphen. This is problematic because it is interpreted as a minus() operation to remove 'data' from 'plink'. You get a NullPointerException since 'plink' hasn't been defined. You can test this easily with the following example:

```
Channel
    .of(1..23, 'X', 'Y')
    .set { chr-names }

chr-names.view()

```

Results:

```
Cannot invoke method minus() on null object

```

The solution is to just use an underscore instead of a hyphen:

```
Channel
    .of(1..23, 'X', 'Y')
    .set { chr_names }

chr_names.view()

```

---

**Follow-up issues:**

I think if you want to use the bfiles plink option (and I think this is reasonable), I think it makes the most sense to just keep your bfiles in a tuple:

```
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bfiles) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bfile "${sample_name}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

```

If you need to reference each bfile separately, you could just unpack the tuple and reference each file accordingly:

```
Channel
    .fromFilePairs( "${params.input_dir}/*.{bed,bim,fam}", size: 3 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .set { plink_data }

process test1{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from plink_data

    output:
    tuple val("${sample_name}-t1"), path("${sample_name}-t1.{bed,bim,fam}") into test1_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.05 \\
        --make-bed \\
        --out "${sample_name}-t1"
    """
}

process test2{

    publishDir params.publish_dir

    module 'plink'

    input:
    tuple val(sample_name), path(bed), path(bim), path(fam) from test1_results

    output:
    tuple val("${sample_name}-t2"), path("${sample_name}-t2.{bed,bim,fam}") into test2_results

    """
    plink \\
        --bed "${bed}" \\
        --bim "${bim}" \\
        --fam "${fam}" \\
        --maf 0.01 \\
        --make-bed \\
        --out "${sample_name}-t2"
    """
}

```

Note that in the first example, `bfiles` could be replaced with `"*"` since no reference is made to the tuple of bfiles anyway. Also note that the list of files produced using the fromFilesPairs operator are returned in lexicographic order.
> 2 votes

---
Tags: nextflow
---


title: differential analysis of chip-seq data

text:
# Question
Title: differential analysis of chip-seq data
I have several sets of chip-seq data. I called the peaks using Macs2. I am pretty new to the field and I will appreciate any help. I wanted to annotate the peaks and see which peaks are shared between them and which ones are different. What is the best way to do this? which package should I use?

Many thanks

# Answer
1. For sharing, you can test the proportion of overlapping peaks using `bedtools intersect` or the `find.overlap` function in `GenomicRanges` packages in R.
2. What makes more sense is perhaps to merge all the peaks together using `bedtools merge`, and then count the number of reads/fragments in each peak for each sample using `featureCounts`. The output of this is a matrix, which is very similar to gene expression results. Then you can do some differential tests on this data.
> 1 votes

# Answer
In a similar question, I recommended that directly comparing called peaks can be somewhat misleading. This topic is addressed in the documentation of the `DiffBind` R package on Bioconductor, in the "Comparison of occupancy and affinity based analyses" section.

But if you're just looking for genome annotation software, the `annotatr` package on Bioconductor is pretty useful
> 0 votes

# Answer
It depends on how many samples you have. If you only compare one sample to the second, you can pool and call peaks and then directly calculate RPKM ratios for each peak. It you have more samples, you may want to use IDR to score differential peaks. You may also use packages such as Fseq or HOMER to treat one group as ChIP experiment while the other group as control to call peaks. Alternatively, you could pool all the samples to call peaks, then calculate read counts in each peak for each sample, and then use DESeq2 or other differential tools to calculate significantly different peaks’ significance and effect size.
> 0 votes

---
Tags: sequence-annotation, sequence-analysis, chip-seq, peak-calling, macs2
---


# Question
Title: Calculate genome coverage and depth from alignment
I have a .bam alignment file and a genome reference .fasta file. I am looking for a easy to use tool (that I can reference in a publication) to calculate the percentage coverage of the reference by the alignment and if possible also: average depth of reads that cover the gene, min/max depths. Running plot would be nice but not essential.

I wrote some python code to calculate coverage percent, but I want to reference the tool I use in a paper and dont want to mess with making a git repo and linking to it in the paper just for this simple analysis.

# Answer
> 2 votes
You could use `samtools coverage` as explained in the manual of samtoools.

Here is a example which is also described on the manual site.

```
samtools coverage -r chr1:1M-12M input.bam

#rname  startpos  endpos    numreads  covbases  coverage  meandepth  meanbaseq  meanmapq
chr1    1000000   12000000  528695    1069995   9.72723   3.50281    34.4       55.8

```

# Answer
> 0 votes
I've used https://genome.sph.umich.edu/wiki/BamUtil before. I think this has all the functionality you want and the Devs are responsive which is nice.

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


# Question
Title: Aggregate sequencing/mapping/etc. metrics from cellranger across Illumina samples
I have a number of single-cell projects processed with cellranger from 10x Genomics. The pipeline produces a number of handy metrics that are summarized for each Illumina sample in web page. These include things like the estimated number of cells, the fraction of reads in cells, and mapping rates. What I would like is to aggregate this information from multiple Illumina samples into a single file. MultiQC is the obvious choice for such a task, but it appears that no one has written a module to process cellranger data for it. Are there any other prepackaged ways to perform this aggregation?

If there's no convenient tool for this I'll just write a MultiQC module, but I'd obviously like to avoid the extra effort unless it's really needed.

# Answer
> 1 votes
I wrote a simple script to do it. Not sure whether that's what you want.

#!/bin/bash
#CellrangerSummary.sh Sample1 Sample2
csvtool transpose $1/outs/*summary.csv | csvcut -c 1 > summary.csv
for sample in "$@"
    do
        csvtool transpose $sample/outs/*summary.csv > tempfilesforcsvtool.$sample
        csvjoin -c 1 summary.csv tempfilesforcsvtool.$sample > temp2
    mv temp2 summary.csv
rm tempfilesforcsvtool.$sample
    done

csvformat -T summary.csv > summary.tsv
rm summary.csv
```

---
Tags: scrnaseq, cellranger
---


title: Download multiple SRA files

text:
# Question
Title: Download multiple SRA files
I want to download all SRA file from the following project. Is there a method to download all the SRA files at the same time?

# Answer
A quick look at your link tells me the SRR numbers run from SRR837819 to SRR837856. You can use fastq-dump from the sratoolkit, and make a for loop around it in bash.

Something like this should work:

```
for (( i = 19; i <= 56; i++ ))
  do
  fastq-dump --accession SRR8378$i
done

```

After reading Devon Ryan's answer, I realize that you asked for SRA files instead of fastq. This can also be done with prefetch of the sratoolkit:

```
for (( i = 19; i <= 56; i++ ))
  do
  prefetch SRR8378$i
done

```
> 4 votes

# Answer
Assuming you ultimately just want the fastq files and you know the SRR (run) numbers, I would download them from here: ftp://ftp.sra.ebi.ac.uk/vol1/fastq/

As for downloading multiple files, I've just used multiple wget commands. I don't know of a way to download all of the files together in like a zipped folder or anything :/
> 3 votes

# Answer
I suggest you follow the advice in Eric A Brenner's answer and just download the fastq files. However, if you really really want to use the SRA files for some reason, note that you can use parallel-fastq-dump to make things faster. Do follow its advice regarding using `prefetch`.

You'd need to combine that with the answer from b.nota (i.e., put the commands in a `for` loop).
> 3 votes

# Answer
I was able to find a solution for this using Entrez Direct and SRA toolkit :) If you have the project number or SRA project number, yours would be SRP022054 in this case for the 36 SRAs, you can use esearch to make a query like so and pipe it into SRA toolkit with this one liner:

`esearch -db sra -query SRP022054 | efetch --format runinfo | cut -d ',' -f 1 | grep SRR | head -5 | xargs fastq-dump --skip-technical --readids --read-filter pass --dumpbase --split-3`
> 1 votes

# Answer
You can download an entire project using pysradb:

```
pysradb download -p SRP022054

```

It retains the same schema as SRA.
> 0 votes

# Answer
Once you install SRA tools, you can grab a BioProject/SRP accession list (Send to - File - Accession List) and save it to default SRA tools /sra directory, then run (with /sra as cwd):

```
prefetch --option-file SraAccList.txt
cat SraAccList.txt | xargs fasterq-dump --outdir "/your-directory/for-fastq"

```

P.S. The text file can be any list of accessions, separated by return. Tested in Bash
> 0 votes

# Answer
This online tool creates download links for the less tech savy:

https://sra-explorer.info/#

Using PRJNA201245 as ID, it provides several codes including a bash script for downloading all files.
> 0 votes

---
Tags: rna-seq, data-download
---


# Question
Title: Nextflow: No signature of Method: chr() is applicable for argument types; unknown variable
I was trying to use plink in nextflow. My first process will split the input (single bed|bim|fam file) into 22 chromosomes (each with bed|bim|fam files) which would be used in subsequent processes. My code and error log is below. Any help?

```
params.dir='data/'
params.publish="output2/"
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}

Channel
.of(1..22)
.set {chr_names}

process step1a {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from data
output:
tuple val("${prefix}"), path("${prefix}.step1a.{bed,bim,fam}") into step1a_results
tuple val("${prefix}"), path("${prefix}.step1b.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.step1b.bim") into step1a_bims
path("${prefix}.log") into step1a_logs
path("${prefix}.step1a.dupvar") into step1a_duplicates
path("${prefix}.step1a.log") into step1a_log
path("${prefix}.step1b.log") into step1_log
script:
prefix=chr"${chr_names}"
"""
plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
"""
}
process step1b {
publishDir params.publish
input:
tuple val(sample_name), path(bfiles) from step1a_out
tuple val(bims), path("${prefix}.step1b.bim") from step1a_bims
output:
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
path("${sample_name}.step1c.log") into step1b_logs
path("${sample_name}.step1b.snplist.txt") into step1b_snplist
script:
"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
plink --bfile "${prefix}.step1b --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""
}

```

The error message is:

```
N E X T F L O W  ~  version 20.10.0
Launching `t4.hf` [zen_picasso] - revision: b0ff8db413
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
[-        ] process > step1a -
[-        ] process > step1b -
Error executing process > 'step1a (1)'

Caused by:
  No signature of method: chr() is applicable for argument types: (org.codehaus.groovy.runtime.GStringImpl) values: [DataflowQueue(queue=[DataflowVariable(value=1), DataflowVariable(value=2), DataflowVariable(value=3), DataflowVariable(value=4), DataflowVariable(value=5), DataflowVariable(value=6), DataflowVariable(value=7), DataflowVariable(value=8), DataflowVariable(value=9), DataflowVariable(value=10), DataflowVariable(value=11), DataflowVariable(value=12), DataflowVariable(value=13), DataflowVariable(value=14), DataflowVariable(value=15), DataflowVariable(value=16), DataflowVariable(value=17), DataflowVariable(value=18), DataflowVariable(value=19), DataflowVariable(value=20), DataflowVariable(value=21), DataflowVariable(value=22), DataflowVariable(value=groovyx.gpars.dataflow.operator.PoisonPill@28c22659)])]
Possible solutions: run(), run(), any(), is(java.lang.Object), use([Ljava.lang.Object;), tap(groovy.lang.Closure)

Source block:
  prefix=chr"${chr_names}"
  """
  plink --bfile "${sample_name}" --chr "${chr_names}" --out "${prefix}" --make-bed
  plink --bfile "${prefix}" --list-duplicate-vars ids-only suppress-first --out "${prefix}.step1a"
  plink --bfile "${prefix}" --exclude "${prefix}.step1a.dupvar" --make-bed --out "${prefix}.step1a"
  plink --bfile "${prefix}.step1a" --snps-only just-acgt --make-bed --out "${prefix}.step1b"
  """

Tip: when you have fixed the problem you can continue the execution adding the option `-resume` to the run command line

```

If I use "prefix=chr"${chr\_names}" " inside the code block (""") it shows: unknown variable 'prefix'

# Answer
The documation re conditional scripts is useful here:

> Process scripts can contain conditional statements by simply prefixing the script block with the keyword `script:`. By doing that the interpreter will evaluate all the following statements as a code block that must return the script string to be executed.

But the statements don't necessarily need to be conditional. I.e. they can just be declarations like in your example. You get the error 'No signature of Method: chr()' because Nextflow interprets 'chr' as a method call when you use:

```
prefix = chr"${chr_name}"

```

If your intention is to simply add the 'chr' prefix to each chromosome name, you'd need:

```
prefix = "chr${chr_name}"

```

Notice that I changed your variable name. This is because in your code, you still need to feed in the channel of chromosome names. The easiest way to do this is to use the each qualifier:

```
process myprocess {

    input:
    tuple val(sample_name), path(bfiles) from plink_data
    each chr_name from chr_names

    script:
    prefix = "chr${chr_name}"

    """
    echo "My sample name is: ${sample_name}"
    echo "My prefixed chromosome name is: ${prefix}"
    """
}

```

Unless you need chromosome names without the 'chr' prefix, another solution might be to simply add the prefix to each one initially and simply avoid having to prefix them later:

```
Channel
    .of(1..22)
    .map { "chr${it}" }
    .set { chr_names }

```
> 4 votes

---
Tags: plink, nextflow
---


# Question
Title: Nextflow: ordering the execution of processes; outputs of process A are inputs of process B
I need to build a pipeline where 2nd process starts after completion of the 1st process and the 3rd after 2nd. As I understand nextflow's processes run parallel. But I need it to run in a order 1st\>2nd\>3rd processes and so on. Because each process runs on an array and produce multiple outputs which are used by the next. If they run paralleled, the 2nd process may look for input 20 whereas the 1st process produced output 10 and then the process stops, because the command in the 2nd process does not find the input. Is there any way to accomplish this. I tried as follows:

```
process A{
input:
tuple val(some_name), path(some path) from some_channel
output:
tuple val(some val), path(some path) into A_some_ch
tuple va(some other va), path(some other path) into A_some_other_ch
val true into done_ch
script:
"""
some command out some files
some other command out some other files
"""
}
process B {
input:
val flag from done_ch
val (some val), path(some path) from A_some_ch
val (some other val), path(some other path) from A_some_other_ch
output:
some val and path
script:
"""
some commands out some files
"""
}

```

I got the idea using of flag from here: https://github.com/nextflow-io/patterns/blob/master/docs/mock-dependency.adoc

But if I use the flag, I got:

```
Operation 'set' does not allow more than one target name

```

Does it clarify the issue? Any Help? Best, Zillur

# Answer
Nextflow uses *channels* (asynchronous FIFO queues) to communicate between processes. The interaction between these processes, and ultimately the pipeline execution flow itself, is implicitly defined by these input and output declarations.

The issue with your code, I think, is the use of multiple input channels. Consider the following example that uses only a single *queue* channel. To avoid declaring the inputs again in process A's output declaration, we can output a key to join on as the first element in the output tuple. In the case below, this is simply the 'some\_name' variable. We can then join up the two channels using the join operator as follows:

Channel
    .fromFilePairs( './data/*.txt', size:1 )
    .map { group_key, files -> tuple( group_key, *files ) }
    .into { input_ch_A; input_ch_B }

process A {

    input:
    tuple val(some_name), path(some_path) from input_ch_A

    output:
    tuple val(some_name), val("procA_${some_name}"), path("procA_${some_name}") \
        into procA_results

    """
    touch "procA_${some_name}"
    """
}

process B {

    input:
    tuple val(some_name), path(some_path), val(some_other_val), path(some_other_path) \
        from input_ch_B.join(procA_results)

    """
    true
    """
}
```

> 1 votes

---
Tags: nextflow
---


# Question
Title: GATK 4.2 using IndexFeatureFile with multiple files
I am trying to use GATK IndexFeatureFile to create ".vcf.gz.tbi" file. But I have 480 files.

For the pilot study I used:

```
!./gatk IndexFeatureFile -I HG00096.chrom22.bam.10061.vcf.gz -O HG00096.chrom22.bam.10061.vcf.gz.tbi

```

Doing this procedure for 480 files is not seem very practical. Is there any way to use all file as an input and it takes all files and gives ".tbi" files as an output?

I thought

```
!./gatk IndexFeatureFile -I *.vcf.gz -O *.vcf.gz.tbi

```

or maybe not using output argument (just input). I did not worked for me.

# Answer
> 2 votes
You could run loop over the files in your shell, running a separate command for each file. For example, in Bash (or another Bourne shell):

```
for file in *.vcf.gz; do ./gatk IndexFeatureFile -I ${file} -O ${file}.tbi; done

```

---
Tags: vcf, gatk
---


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

# Answer
> -1 votes
https://github.com/BrooksLabUCSC/flair was excellent. Highly recommended. One could also ask for TPM.

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

# Answer
> 0 votes
Stringtie now does long reads too. We've found flair to be quite difficult to install and to crash often. That said, stringtie isn't very good about showing you its working (e.g. I have no idea how to get it to give the assembly it uses to figure out all the isoforms). Its values for TPM/FPKM also seem a little weird, in that they're perfectly correlated with each other, and for transcripts they're perfectly correlated with coverage too.

https://genomebiology.biomedcentral.com/articles/10.1186/s13059-019-1910-1 https://github.com/gpertea/stringtie

TALON is another tool that does isoform assembly. I'm not sure whether it does quantification though. https://github.com/mortazavilab/TALON#

---
Tags: rna-seq, nanopore, long-reads
---


# Question
Title: RNASeq analysis using featureCount and EdgeR
I am using a pipeline (bam -\> featurecount-\> EdgeR) to do some RNASeq analysis of several groups and sub-groups. For example, I have the following dataset with two types (T1 and T2) and T1 has three groups G1, G2 and G3, T2 has G5 and G6 groups. I would like to find TPM numbers, fold-change and FDR for all the genes and Groups.

```
Type Group SRA_Run#
T1   G1.   <some SRA ids e.g. SRR123456, SRR112233>
T1   G2.   <some SRA>
T1.  G3.   <some SRA>
T2.  G5.   <some SRA>
T2.  G6.   <some SRA>

```

I have individual BAM files generated for each SRA Run. I also individual read counts from featureCount for each SRA Run. I was not sure if I need to merge individual read counts for each group or run featureCount using all bam files of a group at once. Does it make any difference in the analysis that follows.

I also would like to know if I can run EdgeR at once for all the types and groups at once for doing differential gene expression for each group of Type 1 (T1) against every other group of Type 2 (T2). I am bit stuck at this point how to proceed. Any help would be appreciated.

# Answer
It makes no difference if you process the BAM files one at a time with featureCounts or all together, except that it changes how you have to read the files into R.

You can supply edgeR with lists of contrasts to have it compute fold-changes and p-values for. Please have a look at the edgeR user guide for examples.
> 1 votes

---
Tags: rna-seq, edger, featurecounts
---


title: nextflow: Filter outputs of a process

text:
# Question
Title: nextflow: Filter outputs of a process
How could I filter outputs of a process in the input of the next process? Filtering works fine in channel, but if I try to filter outputs I got compilation error. I tried as follows:

```
Channel
.fromFilePairs("${params.dir}*.{bed,bim,fam}",size:3)
.set {data}
    process step1a {
    input:
    tuple val(sample_name), path(bfiles) from data
    output:
    tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)
    script:
    """
    plink --bfile "${sample_name}" --chr 1-22 --out "${prefix}.step1b" --make-bed
    """
    }
    process step1b {
    input:
    tuple val(bim), path(bims) from step1a_bim.filter{bim, files -> name =~/*bim/}
    tuple val(sample_name), path(bfiles) from step1a_out
    output:
    path("${prefix}.step1b.snplist.txt")
    tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into step1b_results
    script:
    """
    awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${prefix}.step1b.bim" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"
"""

```

It gives me

```
Compilation error
-cause: Unexpected input: '{'
process step1b {

```

Step1a gives me chr1.step1b.bed|bim|fam, chr2.step1b.bed|bim|fam.......chr22.step1b.bed|bim|fam (22 plink files) I need to access prefix of those outputs to run plink and other tools to generate chr1.step1c.bed|bim|fam ...chr22.step1c.bed|bim|fam (another 22 plink files). run awk to generate chr1.step1b.snplist.txt .... chr22.step1b.snplist.txt (22 txt files). Any help to filter the inputs so that I could use only \*bim files from step1a\_bim? Also how to select prefix of files in the input of processes that are from step1a\_out to run plink? Does it make sense now?

Best Reagrds Zillur

# Answer
Instead of outputting all files and filtering for only bim, output only bims:

```
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into (step1a_out,step1a_bim)

```

Do

```
output:
tuple val("${prefix}"), path("${prefix}.{bed,bim,fam}") into step1a_out
tuple val("${prefix}"), path("${prefix}.bim") into step1a_bim

```

So you can skip your complicated filter:

process step1b {
    input:
    tuple val(bim), path(bimfile) from step1a_bim
    output:
    path("${prefix}.step1b.snplist.txt")
tuple val("${sample_name}"), path("${sample_name}.step1c.{bed,bim,fam}") into
step1b_results
    script:
    """
    awk '{ if ((\\\\\$5=="T" && \\\\\$6=="A")||(\\\\\$5=="A" && \\\\\$6=="T")||(\\\\\$5=="C" && \\\\\$6=="G")||(\\\\\$5=="G" && \\\\\$6=="C")) print \\\\\$2, "ambig" ; else print \\\\\$2 ;}' "$bimfile" | grep ambig > "${prefix}.step1b.snplist.txt"
    plink --bfile "${prefix}.step1b" --exclude "${prefix}.step1b.snplist.txt" --make-bed --out "${prefix}.step1c"

"""

```

> 1 votes

---
Tags: plink, nextflow
---


title: Un-normalize DESeq2 counts

text:
# Question
Title: Un-normalize DESeq2 counts
I obtained a matrix of RNA-seq count data that has been normalized by DESeq2's median of ratio method.

I know that DESeq2 wants to take in un-normalized counts, but I do not have access to those data.

How do I best proceed here if I want to perform DEG analysis using DESeq2?

I know I can always start from .fastq files, but that would be so much extra work.

I don't think I can un-normalize the data because I don't know the scaling factors.

Any suggestions?

# Answer
You can read in the normalized count table and don't normalize the data, but my advice here is not to do that. The rounding of the normalized matrix introduces some noise, but I think the larger issue is how sure are you that the table you are working with, is exactly a count table of normalized counts from DESeq2 ?

It really doesn't take much time to align and get count table again, rather than going back and forth wondering whether you did the normalization correctly.

I can show below with an example how it will look like, first we make an example dataset and obtain the normalized counts:

```
library(DESeq2)
set.seed(111)
sz = runif(6,min=0.5,max=1.5)
x = makeExampleDESeqDataSet(sizeFactors=sz,m=6)
x = estimateSizeFactors(x)
ncounts = counts(x,normalize=TRUE)
x = estimateDispersions(x)

```

Now suppose we feed this into a new DESeq2 object, but before that you need to do is check most of the size factors are close to 1:

```
estimateSizeFactorsForMatrix(round(ncounts))
 sample1  sample2  sample3  sample4  sample5  sample6
1.039179 1.045991 1.048276 1.045872 1.032420 1.052006

```

If you see something above that suggest a deviation from 1, something is wrong from the normalization and do not even think about using it.

Below you can use the normalized counts and you can see the difference is not huge:

```
dds = DESeqDataSetFromMatrix(round(ncounts),colData(x),~condition)
sizeFactors(dds) = 1
dds = estimateDispersions(dds)

```

Compare the new count table:

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

```

And your input:

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

```

We can check the dispersion, for lowly expressed genes you see some difference, but on the whole difference

```
with(elementMetadata(x),plot(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#FFC93C80"))
with(elementMetadata(dds),points(log10(baseMean),dispGeneEst,pch=20,cex=0.7,col="#DBF6E980"))

```

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

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

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


title: Telling grep to treat N as [ATCG]

text:
# Question
Title: Telling grep to treat N as [ATCG]
Okay so I'm using grep to try and get a preview of some trimming operations that are not going as expected.. Lets say that my sequence in the FastQ file is: ATNGCNATCG

What I want to do is..

`grep "ATCGCTATCG" my.fastq`

..and match the sequence given above

Surely there is some way or some existing tool that I can use besides doing: `grep "[A|N][T|N][C|N][G|N][C|N]...etc."`

# Answer
If you want to stick to `grep`, use a scripting language such as Perl to generate the regex programmatically. For example:

```
perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' ATCGCTATCG

```

Prints:

```
[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]

```

You can use it in `grep` like so:

```
grep '[AN][TN][CN][GN][CN][TN][AN][TN][CN][GN]' <<< ATNGCNATCG

```

Prints:

```
ATNGCNATCG

```

If that works for you, you could make it into a little bash function that also runs the `grep`. Add these lines to your `~/.bashrc`:

```
grepN(){
    seq="$1"
    file="$2"
    pattern=$(perl -le 'print join "", map "[${_}N]", split //, $ARGV[0];' "$seq")
    grep "$pattern" "$file"
}

```

You can now run

```
grepN ATCGCTATCG my.fastq

```

Of course, this is not a good idea since the sequence might be in different lines, but that's what you were doing originally.
> 5 votes

---
Tags: fastq, trimming, grep
---


# Question
Title: How to find genes whose expression vary the most in time-series?
Hi I am reading this paper and trying to analyse the same datasets in my work: **dynGENIE3: dynamical GENIE3 for the inference of gene networks from time series expression data**

In this paper, the authors are working on several time course gene expression datasets, including "Drosophila melanogaster embryogenesis" dataset (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE6186) . This dataset has the whole genome ( 14064 genes) of Drosophila melanogaster measured over 31 time-points.

In the dynGenie3 paper, we analyse a subset of these 14064 genes, as they write: *"We focused our analysis on the 1000 genes whose expression vary the most across the time series."*

I want to know how do we measure the variance in expression across time series? What statistical tests might have been performed? I couldn't find any more information in the paper so I am assuming it might be something simple and frequently done. Need some hints.

# Answer
> 0 votes
Maybe my thoughts are too simple, but wouldn't the computation of the variance per gene suffice for this? You really do not need a statistical test, since they are not talking about DEGs, just high-variant genes.

---
Tags: gene, gene-expression
---


# Question
Title: How to know the number of SNP in my samples after running the population module in stacks
I am new to stacks for analyzing RAD-Seq data, I have run all the modules. I will like to know to get the total SNP from my sample after running the population module. Is there a command line to use and inside which file?

command:

```
populations -P -M (denovo output) (popmap) -t 5 --vcf --fasta-samples --fasta-loci -O (pop-output) -r 0.8 -p 5

```

Thanks

# Answer
> 1 votes
**Update**

Based on the newly provided command (which has `--vcf`), there should be a VCF file with suffix `.vcf` in the output directory.

**Original answer**

It's pretty hard to answer without more information. For example, it's helpful to link to the software and the module that you are working from (I had never heard of Stacks). What do you mean by "total SNP"? Also, what is the actual command that you ran? There are a lot of possible options!

But (based only on reading the docs for the command), it looks like there are a variety of formats in which the SNP data can be output:

```
File output options:
--ordered-export — if data is reference aligned, exports will be ordered; only a single representative of each overlapping site.
--fasta-loci — output locus consensus sequences in FASTA format..
--fasta-samples — output the sequences of the two haplotypes of each (diploid) sample, for each locus, in FASTA format.
--vcf — output SNPs and haplotypes in Variant Call Format (VCF).
--genepop — output results in GenePop format.
--structure — output results in Structure format.
--radpainter — output results in fineRADstructure/RADpainter format.
--plink — output genotypes in PLINK format.
--hzar — output genotypes in Hybrid Zone Analysis using R (HZAR) format.
--phylip — output nucleotides that are fixed-within, and variant among populations in Phylip format for phylogenetic tree construction.
--phylip-var — include variable sites in the phylip output encoded using IUPAC notation.
--phylip-var-all — include all sequence as well as variable sites in the phylip output encoded using IUPAC notation.
--treemix — output SNPs in a format useable for the TreeMix program (Pickrell and Pritchard).

```

I would suggest examining the command that you used, specifically what you passed as the output format. Here are the example commands from the docs:

```
populations -P dir -b batch_id [-O dir] [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)
populations -V vcf -O dir [-M popmap] (filters) [--fstats] [-k [--window_size=150000] [--bootstrap [-N 100]]] (output formats)

```

Can you track down what you specified as output? That will contain the SNP data / genotypes.

**Possibly you will need to rerun the command while specifying a different, preferred format.**

---
Tags: snp, phylogenetics
---


# Question
Title: Demultiplexing FASTQ file without index information
I am trying to understand how data that was uploaded to SRA (https://www.ncbi.nlm.nih.gov/sra?LinkName=biosample\_sra&from\_uid=4510743) can be analyzed with the assumption that the FASTQ file should contain replicate information (3 replicates within the sample). I downloaded the file and ran it with Nextflow, and it was mapped to the human genome and determined to contain unstranded reads of 50bp, both directions.

Which steps can I take to infer if the FASTQ file indeed is not demultiplexed?

Many thanks, Geo

The header of the FASTQ looks like this:

```
@SRR3467208.1 1/1
NTCTCTATGTCCACTCTGGAGCCTTTAAGTGCCACTTGAGGGCCTTTAAC
+
#4:D;DDEHHFHHIJJIJJGHGHJJI@HFIIJJHJJIG?9CFHIIIHHHG
@SRR3467208.2 2/1
NGTGCTGTGTGTGCATGTGTGTGCGTGTGTGTGCTGTGCGTTTGTGTGTG
+
#1:B:BDDFFFFFFFBGFHIFIIFGFFGFGFGIIIIIIEFGGIIFFIFII
@SRR3467208.3 3/1
NGTCACTTGGCATCTGGCCATCGGGCTGGATGCCGTGTTCCAGGCAGTAG
+
#1:A=DDDB<<CFAFGI@FGG9A@68CEG;DB?<DGGGE49?BBF2B=FF
@SRR3467208.4 4/1
NTCCGGGGCTCCAGCAACCAGAAGAGGGAAAAGCTGTCTTCGGTGATGTC
+

```

# Answer
> 1 votes
I would suggest looking at sample metadata:

> Library:

> Name: TruSeq RNA

> Instrument: Illumina HiSeq 2000

> Strategy: RNA-Seq

> Source: TRANSCRIPTOMIC

> Selection: PolyA

> Layout: PAIRED

You can see that it's paired. You can look at the spot descriptor at the page you link to see more information about the design.

There is no multiplexing information available there. I would instead look at the overall project and its metadata to find the different replicates. It looks like there are 2 samples x 3 replicates each, deposited independently. So it is likely to already be demultiplexed, and you can just download them individually.

"Unstranded" isn't really a meaningful term for this, if I understand correctly. Forward and reverse read pairs don't correspond to forward and reverse on the chromosome, but to the way that sequencing-by-synthesis happens on the machine. So the mapping orientation is not really relevant.

I infer that you are working from a single FASTQ file. If you are working from a single FASTQ file for paired-end out of SRA, then the FASTQ file is not yet resolved into F+R reads. I would suggest rerunning the FASTQ conversion with `--split-3` or similar option.

---
Tags: rna-seq, fastq, demultiplexing
---


# Question
Title: Challenging benchmarks for supervised learning on sparse scRNA-seq data
One challenging aspect of modeling scRNA-seq data is data sparsity, that is, scRNA-seq measurements typically suffer from large fractions of observed zeros (i.e. *dropouts*), where a given gene in a given cell has no unique molecular identifiers or reads mapping to it \[\].

What would be an appropriate dataset to benchmark supervised learning algorithms that operate on single-cell RNA-seq data? In particular, I am looking for challenging datasets (preferably with labels, i.e. cell type classification) where existing methods struggle due to the sparsity of the data. Additionally, does it ever happen that the *dropouts* are confounders of the labels? (i.e. for a given cell type A the *dropout* probability of a gene X is higher than for another cell type B).

Any pointers are welcome.

# Answer
> 2 votes
You could try the PBMC 3k data from the Satija lab:

https://satijalab.org/seurat/articles/pbmc3k\_tutorial.html

For well-annotated data, there's the Single-cell proteo-genomic reference map:

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

Data for that are available here.

If these datasets are too large (i.e. too good), you can use the existing counts as weighting for a random removal of reads. Genes with high expression are likely to be present regardless of sequencing depth, hence the count-adjusted removal.

FWIW, I don't think that supervised learning is a good fit for dealing with data sparsity in cDNA read counts. There are plenty of existing attempts to solve the issues by adjusting the assumed models, including research by the authors of commonly-used bulk differential expression tools (e.g. here).

---
Tags: rna-seq, scrnaseq, single-cell, machine-learning, benchmarking
---


# Question
Title: Nextflow: Runs without error, produce few but not all output files into output channel
I was trying to develop a pipeline using nextflow and plink. My code runs without any error, but got only two output files instead of 22 (one for each chromosome). Here is part my code block, any help?

```
process step5 {
publishDir params.publish, mode: 'copy'
input:
val(hrc) from HRC_1000G_check_bim
val(hrc_site) from hrc_sites_file
tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
output:
tuple val("${prefix}"), val(chr_name), path("LOG-${prefix}${chr_name}.step4-HRC.txt") into step4_log
tuple val("${prefix}"), val(chr_name), path("Exclude-${prefix}${chr_name}.step4-HRC.txt") into step4_exclude
tuple val("${prefix}"), val(chr_name), path("Chromosome-${prefix}${chr_name}.step4-HRC.txt") into step4_chromosome
tuple val("${prefix}"), val(chr_name), path("Position-${prefix}${chr_name}.step4-HRC.txt") into step4_position
tuple val("${prefix}"), val(chr_name), path("Strand-Flip-${prefix}${chr_name}.step4-HRC.txt") into step4_strand
tuple val("${prefix}"), val(chr_name), path("Force-Allele1-${prefix}${chr_name}.step4-HRC.txt") into step4_allele
tuple val("${prefix}"), val(chr_name), path("ID-${prefix}${chr_name}.step4-HRC.txt") into step4_id
script:
"""
perl "${hrc}" -b "${prefix}${chr_name}.step4.bim" -f "${prefix}${chr_name}.step4.frq" -r "${hrc_site}" -h
"""
}

```

The output files are:

```
gsa-out/Chromosome-chr22.step4-HRC.txt     gsa-out/Force-Allele1-chr3.step4-HRC.txt  gsa-out/Position-chr22.step4-HRC.txt
gsa-out/Chromosome-chr3.step4-HRC.txt      gsa-out/ID-chr22.step4-HRC.txt            gsa-out/Position-chr3.step4-HRC.txt
gsa-out/Exclude-chr22.step4-HRC.txt        gsa-out/ID-chr3.step4-HRC.txt             gsa-out/Strand-Flip-chr22.step4-HRC.txt
gsa-out/Exclude-chr3.step4-HRC.txt         gsa-out/LOG-chr22.step4-HRC.txt           gsa-out/Strand-Flip-chr3.step4-HRC.txt
gsa-out/Force-Allele1-chr22.step4-HRC.txt  gsa-out/LOG-chr3.step4-HRC.txt

```

I am not sure why I am getting only 2, instead of 22 output files. With similar output channel previously generated 22 output files in the same script. For eaxmple:

```
process step1b {
publishDir params.publish, mode: 'copy'
input:
tuple val(prefix), val(chr_name), path(bfiles) from step1a_out
output:
tuple val("${prefix}"), val(chr_name), path("${prefix}${chr_name}.step1c.{bed,bim,fam}") into step1b_results
path("${prefix}${chr_name}.step1c.log") into step1b_logs
path("${prefix}${chr_name}.step1b.snplist.txt") into step1b_snplist
script:

"""
awk '{ if ((\\\$5=="T" && \\\$6=="A")||(\\\$5=="A" && \\\$6=="T")||(\\\$5=="C" && \\\$6=="G")||(\\\$5=="G" && \\\$6=="C")) print \\\$2, "ambig" ; else print \\\$2 ;}' "${bfiles[1]
}" | grep ambig > "${prefix}${chr_name}.step1b.snplist.txt"
plink1.9 --bfile "${prefix}${chr_name}.step1b" --exclude "${prefix}${chr_name}.step1b.snplist.txt" --make-bed --out "${prefix}${chr_name}.step1c"
"""
}

```

# Answer
I suspect your `HRC_1000G_check_bim` and `hrc_sites_file` channels are queue channels when they should be value channels. Please again see how multiple channels work. Also, third-party scripts, like your Perl script, would be better placed in a directory called `bin/` at the root of your project repository. Nextflow adds this directory to your $PATH, so scripts in this directory will be automatically accessible in your pipeline without needing to localize them using a separate channel. This way also avoids having to specify an absolute path to invoke them. Just make sure your script is executable: `chmod +x your_script.pl`

```
hrc_sites = file( params.hrc_sites_file )

process step5 {

    input:
    tuple val(prefix), val(chr_name), path(bfiles) from step4_frq
    path hrc_sites

    """
    your_script.pl \\
        -b "${prefix}${chr_name}.step4.bim" \\
        -f "${prefix}${chr_name}.step4.frq" \\
        -r "${hrc_sites}" \\
        -h
    """
}

```
> 1 votes

---
Tags: plink, nextflow
---