title: Counting reads within the intron

text:
# Question
Title: Counting reads within the intron
I have paired-end RNA-seq samples and a list of intron coordinates in bed format. I want to count the intronic reads such that:

```
1. The reads should overlap the intron by atleast 25 bp
2. If both the reads of the pair are within the intron, count it as 1.

```

What tool can I use to compute it.

Thanks

# Answer
> 3 votes
You could use BEDOPS `bedmap` to map reads to introns, using 1) the `--count` operator to do counting of reads overlapping by your criteria; and, 2) the `--indicator` operator to do a true/false operation, where reads are contained entirely within the intron.

For instance, to count reads that overlap introns by at least 25 bases, use `--count` to count and `--bp-over 25` to specify the minimum overlap threshold:

```
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

```

To filter for introns where reads are contained entirely within the intron, use `--indicator` to get introns that meet the criteria of full containment of reads, via `--fraction-map 1`:

```
$ bedmap --indicator --echo --fraction-map 1 --delim '\t' introns.bed reads.bed | awk '($1==1)' | cut -f2- > answer.indicator.bed

```

Once you have `answer.counts.bed` and `answer.indicator.bed`, you can compare those two results with `bedmap --indicator --exact` to see where there are commonalities (true/1) and differences (false/0), and increment the final count where there are identical introns, e.g.:

```
$ bedmap --indicator --echo --exact --delim '\t' answer.counts.bed answer.indicator.bed | awk '{ $NF += $1; print $0; }' | cut -f2- > answer.final.bed

```

If your reads are in a different format, i.e. BAM, you can use `bam2bed` in a process substitution, replacing `reads.bed` with `<(bam2bed < reads.bam)`, e.g.:

```
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed <(bam2bed < reads.bam) > answer.counts.bed

```

Etc.

You can also just convert separately, so that you only need to do it once:

```
$ bam2bed < reads.bam > reads.bed
$ bedmap --echo --count --bp-ovr 25 --delim '\t' introns.bed reads.bed > answer.counts.bed

```

Etc.

# Answer
> 1 votes
You can almost do what you want using `featureCounts` from the `subread` package. The only problem is that featureCounts requires GTF rather than bed.

If you had a GTF of your introns, you could do:

```
featureCounts -a introns.gtf my_bam.bam -f -p --minReadOverlap=25 -o counts.tsv

```

---

The problem with this then, is getting the `GTF` when you currently have a `BED`. One way to get the `GTF` would be to use `CGAT`.

```
cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons -S genes.gtf.gz
cgat gtf2gtf -I genes.gtf.gz --method=exons2introns -S introns.gtf

```

or you could put this together into one pipeline:

```
cgat gtf2gtf -I ensembl.gtf.gz --method=merge-exons  | cgat gtf2gtf --method=exons2introns -S introns.gtf

```

This will first merge all the exons of a gene together to make a single, collapsed transcript. It will that invert this to find that parts that are in the introns of this collapsed transcript. What you end up with is parts of the gene that are never exonic in any of the transcripts.

alternatively, you could convert your bed regions into GTF:

```
cgat bed2gff -I introns.bed -S introns.gtf --as-gtf

```

---
Tags: rna-seq
---


# Question
Title: Is there a tool that can take a protein's amino acid sequence and would display it's locus on the genome?
I have the UNIPROT IDs, PDB IDs and FASTA files of several known proteins. I am looking for a tool that can take as input the protein's amino acid sequence and display the coding nucleotides of those amino acids. Moreover, I look for a tool that would be able to display those nucleotides' location on the human genome.

The final outcome I desire to achieve is to find SNPs in the exons of the protein.

A tool that can handle several amino acid sequences at once is preferred for I have plenty of proteins. I've tried to use UCSC's BLAT human genome search but it was too cumbersome. I've also tried BLAST, but I'd like a more "straightforward" method for displaying the genome locus for each protein. By the way, I've tried to use UCSC's Table Browser but couldn't query several genes simultaneously.

Thank you for your assistance.

# Answer
> 2 votes
I think the mistake here is thinking that you need to use the sequence in order to find the genome co-ordinates. If these proteins come from a well annotated speices, you can easily do this with BioMart.

1. Go the the ENSEMBL web site.
2. From the top menu, select BioMart
3. Under "-CHOSE DATABASE-" select Ensembl Genes, and under "-CHOOSE DATASET-" select your species (eg. Human Genes).
4. Select "Filters" from the lefthand menu and then expand the "GENE" section.
5. Select the "Input external references ID list \[Max 500 advised\]" and from the drop-down select "PDB IDs" or the correct "Uniprot KB" ID.
6. Paste your IDs into the box.
7. Select "Attributes" on the left-hand menu. Making sure you've the "Fetaures" radio box select, you can select Chromosome and either Gene Start/Gene End or Transcript start/Transcript end from the GENE section. These coordinates will also contain the untranslated regions unfortunately. You may also want to play about with the various things under "PROTEIN DOMAINS AND FAMILIES", particurlarly under "Protein Features", although I'm not sure what they all mean.
8. When you are happy with the selection, select the big "Results" at the top to get the coords.

You can also get out Sequences for both the full transcript and the CDS and known germline and somatic variants that fall within the boundaries of the genes.

# Answer
> 1 votes
Have you tried BLAST, or specifically tBLASTn? This will take a protein and compare it any set of reference genomes. You can also BLAST multiple sequences at once.

---
Tags: fasta, sequence-alignment, snp, pdb, uniprot
---


# Question
Title: How to demultiplex pair-end fastq reads with barcode 2 in the identifier line?
I have multiplexed pair-end fastq reads with dual barcodes. The issue is that one barcode is present in the header and one is present at the beginning of the read. I need a method to demultiplex this data, but in order to assign a read to an individual, both barcodes are required, as there is overlap between the barcodes. It seems there are packages available to demultiplex using header ID or in-line barcodes to demultiplex, but not both.

example reads:

```
@700819F:525:HT235BCXX:2:1101:1139:2144 1:N:0:ATCACGAT
CGAATTGCAGATTTTTTCTGAATAAAGCAGTGCAATAAAATTCCCCGCAAAAACACTTNANNNGNNNNNNNNNNNNNNNNNNNNANNNNNNGNTAATAAA
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIII#<###.####################.######<#.<<GGGI
@700819F:525:HT235BCXX:2:1101:1212:2172 1:N:0:ATCACGAT
AAGGATGCAGGGCATCTCCCTCAGGCTGCGCTCTATCGAAGTCATCCCAGAATTAGATTCCGACCACAGACCAGTCTTAGTCAAACTAGGACCCGAGTGT
+
GGGGGIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIIGIIIIIIIIIIIIIIIIIIIIGIIIIIIGGIIIGIIIIIIIIIIIIIGGIGIGIIIGIGI
@700819F:525:HT235BCXX:2:1101:1110:2173 1:N:0:ATCACGAT
GGTTGTGCAGAAAGAGTTGCTGATAAACTTAGCCATGCAGAACAGAATTATGAGTTAGAAGTATGTATATATATACCAATCACTATATCAACCCATTACC
+
<G.G<GGIIIG.GA.GAGAAGG<AA<A<.<<<GA<.<<.G<.G<<A<GGAA....<G.G..<<.GA..<A.<GG<<<.<..<<.GG.A..G..<<<.<<G

```

One of the barcodes is the sequence found at the end of the header line (ATCACGAT). The second barcode is found at the beginning of the sequence (CGAAT in read 1). Both are required for assignment to an individual.

# Answer
You could try something like that:

```
#!/bin/sh

fq_in="${1}"

bioawk -c fastx '{print substr($comment, 7, 15)}' ${fq_in} | uniq \
    > "barcode1.txt"
bioawk -c fastx '{print substr($seq, 1, 5)}' ${fq_in} | uniq \
    > "barcode2.txt"

for bc1 in $(cat "barcode1.txt")
do
    for bc2 in $(cat "barcode2.txt")
    do
        fq_out="${bc1}_${bc2}.fq"
        bioawk -c fastx -v bc1=${bc1} -v bc2=${bc2} \
            'substr($comment, 7, 15) == bc1 && substr($seq, 1, 5) == bc2 {print "@"$name" "$comment"\n"$seq"\n+\n"$qual}' \
            ${fq_in} \
            > ${fq_out}
    done
done

```

This might be slow. You can a little bit of time if you directly provide lists of barcodes for the for loops.

It could be more efficient to write a custom program that can write to multiple files. It happens that I wrote something that does this in python just yesterday, but only reading the inline barcode. Modifying it to also use the in-header barcode is likely simple.
> 2 votes

# Answer
Since you've found separate programs to demultiplex by Illumina barcode and in-line barcode (I presume this is actually just the sequence from the cut site), just run one and then the other on its output. Typically the Illumina barcode would have been used to demultiplex the files before you even received them (if not, your sequencing facility is incompetent), in which case you just need to demultiplex by the in-line barcode
> 1 votes

# Answer
The program `extract` from `UMI-Tools` will extract the barcode from the read and add it to the ID. It will also handle pair-end data so that even if the barcode is only on one of the reads, it will be added to the ID of both reads.
> 1 votes

---
Tags: fastq
---


# Question
Title: How to concatenate "by chromosome"-VCFs?
I have a several VCFs which are VCF which only contain information by chromosome. That is, there's a chromosome 1 VCF (with only chr1), a chromosome 2 VCF (with only chr2), etc.

I checked to make sure that these VCFs were valid via `VCFtools`, i.e.

```
$ vcf-validator chr1.vcf

```

which works---these are valid VCFs I was given.

Now, I would like to combine these VCFs into one VCF.

I naïvely tried the following `cat` operation:

```
$ cat chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf > total_chroms.vcf

```

This doesn't work properly though

```
$ vcf-validator total_chroms.vcf
The header tag 'contig' not present for CHROM=chr1. (Not required but highly recommended.)
Could not parse the line, wrong number of columns: [##fileformat=VCFv4.2\n]
 at /path/vcftools-0.1.14/perl/Vcf.pm line 172, <__ANONIO__> line 191016.
     Vcf::throw('Vcf4_2=HASH(0x1ae7208)', 'Could not parse the line, wrong number of columns: [##filefor...') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 335
     Vcf::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 3457
     Vcf4_1::next_data_array('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14/perl/Vcf.pm line 2574
     VcfReader::run_validation('Vcf4_2=HASH(0x1ae7208)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 60
     main::do_validation('HASH(0x16ada68)') called at  /path/vcftools-0.1.14//bin/vcf-validator line 14
$

```

What tools are available to merge these VCF's together into a total VCF?

# Answer
I would recommend bcftools concat. You can't just `cat` them together because each file has a header section. The bcftools command will handle all that for you. Each vcf file must be sorted prior to calling `concat`

`bcftools concat -o total_chroms.vcf chr1.vcf chr2.vcf chr3.vcf ... chrX.vcf`
> 11 votes

# Answer
The best tool for the job is probably `bcftools` as suggested by Bioathlete, but you can also do this manually. You just need to collect all header lines from all vcf files, remove any duplicates, and then print all the headers + the actual data into the new one:

```
grep '^##' chr*vcf | sort | uniq > all.vcf
grep -m1 '^#CHR' chr1.vcf >> all.vcf  ## Get the chr header line
grep -v '^#' chr*vcf >> all.vcf

```
> 2 votes

# Answer
use picard **GatherVcf**: http://broadinstitute.github.io/picard/command-line-overview.html#GatherVcfs

> Gathers multiple VCF files from a scatter operation into a single VCF file. Input files must be supplied in genomic order and must not have events at overlapping positions.
> 2 votes

---
Tags: vcf, vcftools
---


# Question
Title: bed entries collapse based on score column
I have a bed file where some entries have the exact same chr:start-end but change in the name and score column. I would like to apply bedtools or a similar tool to obtain a filtered file where for each set of chr:start-end entries, only the one with the highest score (column 5) is used, the others filtered out. I would like to do this with a bed format toolkit (bedtools, etc) rather than scripting it.

# Answer
> 2 votes
Use BEDOPS `bedmap --max-element` or `bedmap --max-element-rand` to grab the maximal-scoring element over mapped intervals, e.g.:

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

```

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

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

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

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

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

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

# Answer
> 1 votes
groupby from bedtools might help you:

```
bedtools groupby -i variantsToRepeats.bed -grp 1-3 -c 5 -o max

```

`-grp` to select the columns to group by
`-c` to select the columns to summarize
`-o` to select how to summarize the column selected.

---
Tags: bed, bedtools
---


# Question
Title: Make ipyrad use cuda-enabled NVidia card on Ubuntu
I want to use ipyrad on a new Ubuntu machine that has an NVidia Quadro K2000 card with 384 cores. One can configure ipyrad to run on a linux cluster. Do I have any options to get ipyrad to access these cores?

I'm assuming that I will be able to install cuda drivers on the machine. If ipyrad can't be made to do this, is there another way to do de novo assembly of GBS data that will take advantage of this resource?

# Answer
> 2 votes
You would have to get CUDA Python set up on your system and likely modify ipyrad to actually use it (for example, by incorporating numba). So while this is possible, it will require significant effort on your part. In general, there are very few bioinformatics tools that support GPUs.

---
Tags: python, assembly, genotyping, radseq
---


# Question
Title: What is the --rev_comp_mapping_barcode parameter in the QIIME1 script extract_barcodes.py?
I'm using the QIIME1 scripts extract\_barcodes.py to extract barcodes from a dual-indexed MiSeq library, and then demultiplex my reads using split\_libraries\_fastq.py. I have been having trouble correctly using a couple of parameters relating to reverse complementing barcodes.

When extracting the barcodes with extract\_barcodes.py, I use `--rev_comp_bc1` (as directed by the individual in charge of sequencing). Then when I demultiplex the reads with split\_libraries\_fastq.py, sometimes I need to include the parameter `--rev_comp_mapping_barcodes` to get the script to work, and sometimes I don't. I'm not sure why, unless the barcodes I'm receiving in my mapping files are sometimes written in different orientations (I don't write them).

I think the parameter `--rev_comp_bc1` is self-explanatory, but could someone explain what `--rev_comp_mapping_barcodes` does and when to use it?
To make things simple we can assume that the two barcodes are written from 5' to 3' in the mapping file.

# Answer
> 1 votes
As the help page says:

> --rev\_comp\_mapping\_barcodes
> Reverse complement barcode in mapping before lookup (useful if barcodes in mapping file are reverse complements of golay codes) \[default: False\]

As the help page of extract\_barcode says:

> --rev\_comp\_bc1
> Reverse complement barcode 1 before writing \[default: False\]

You seem to be doing two complete different things. extract\_barcode.py is designed to

> format fastq sequence and barcode data so they are compatible with split\_libraries\_fastq.py

While the first one is for:

> demultiplexing of Fastq sequence data where barcodes and sequences are contained in two separate fastq files

So it seems that depending on the machine used or the previous processing you are given with the barcode in another file or together. I would advice to ask the person who gave you this data.

---
Tags: ngs, sequence-analysis, quality-control, barcode
---


# Question
Title: Finding the nearest gene to a specified gene region
I've got some data on genomic regions of interest in a bedgraph-like format from the *Mus musculus* \[mm10\] genome, and would like to find out the nearest gene to these regions.

Some of the regions will reside within a gene (in which case I'd like to know that they're *within* a gene), some regions may reside within multiple genes (in which case I'd like to know all gene names), and some regions are not within a gene (in which case I'd like to know the closest gene: the shortest sequence gap between either end of the gene and either end of the tagged genomic region. In this case I want to do an unstranded search, because it's possible that I might have regions describing suppressing transcripts. Here's some example data:

```
Ref         Start             End      fwd     rev
chrX     73716152        73716152        0      -2
chrX    167207094       167209162        0      -3
chrY     30770844        30772724        3       0
chr1     24613189        24613641        0       3
chr1     24613971        24614812        0       4
chr1     35788729        35789387       -5       0
chr1     93151586        93160225        0      -2
chr1    107511454       107525597       -2       0
chr1    130882067       130887413        0      -3
chr1    134182466       134189942        2       0
chr1    135818954       135832956       -2       0
chr1    149520498       149521157        0      -3
chr1    166166500       166166500       -2       0
chr2      5379059         5379374        0      -2
chr2     32379110        32381837        0      -2
chr2     32381854        32381858        0      -2
chr2     32384633        32384633        0      -2
chr2     32384636        32387773        0      -3
chr2    120539390       120563799        0      -2
chr2    125946785       125975527        2       0
chr2    127633311       127633684        0      -2
chr2    127634841       127656210        0      -3
chr2    164562717       164562717       -2       0
chr2    164562719       164562720       -2       0
chr2    164568498       164568499       -2       0
chr2    181188007       181188388       -2       0

```

I'd like to add another 'Gene' column to this table that gives the official names of the nearest gene(s) for each region defined by `<Ref>:<Start>-<End>`.

# Answer
You're looking for `bedtools closest`, which thankfully defaults both considering overlapping genes and outputting all hits in case of ties.

```
  bedtools closest -a foo.txt -b genes.bed | awk '{foo = sprintf("%s:%i-%i", $6, $7, $8); OFS="\t"; print $1, $2, $3, $4, $5, foo}'

```

You'll need to strip the header and ensure the file is sorted. I've piped the output to `awk` to reformat things. The output then is something like:

```
chr1    24613189    24613641    0   3   chr1:24613188-24613971
chr1    24613971    24614812    0   4   chr1:24613973-24614651
chr1    35788729    35789387    -5  0   chr1:35806973-35810462
chr1    93151586    93160225    0   -2  chr1:93151348-93160948
chr1    93151586    93160225    0   -2  chr1:93151582-93160870
chr1    107511454   107525597   -2  0   chr1:107511422-107525598

```

Those results are for the Gencode M15 release (the BED file I used has transcripts rather than genes, so you'll get a bit different coordinates). As an aside, the gene name would normally be output as well, just add column 9 to the output.
> 4 votes

# Answer
With the prompting from @Devon\_Ryan's answer, I came up with a solution that worked well for me. I downloaded the mm10 RefSeq gene locations from UCSC \[Fields: chrom, txStart, txEnd, name2\], and sorted both input files using `bedtools sort`. From there the `bedtools closest` command got me almost to where I wanted to go:

```
(echo -e "Chr\tStart\tEnd\tfwd\trev\tm_Chr\tm_Start\tm_End\tGene\tDistance";
  bedtools closest -d -a regions.bed -b NCBI_RefSeq_mm10.bed) > annotated.tsv

```

Output example:

```
Chr     Start   End     fwd     rev     m_Chr   m_Start m_End   Gene    Distance
chrX    73716152        73716152        0       -2      chrX    73686177        73716204        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73717815        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716401        Bcap31  0
chrX    73716152        73716152        0       -2      chrX    73686177        73716502        Bcap31  0
chrX    167207094       167209162       0       -3      chrX    167207092       167213479       Tmsb4x  0
chrX    167207094       167209162       0       -3      chrX    167207093       167209218       Tmsb4x  0
chrY    30770844        30772724        3       0       chrY    30725674        30738546        Gm29582 32299
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613189        24613641        0       3       chr1    24257682        24587535        Col19a1 25655
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    24613971        24614812        0       4       chr1    24257682        24587535        Col19a1 26437
chr1    35788729        35789387        -5      0       chr1    35869594        35881164        1110002O04Rik   80208
chr1    93151586        93160225        0       -2      chr1    93151353        93160935        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151353        93160957        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160948        2310007B03Rik   0
chr1    93151586        93160225        0       -2      chr1    93151354        93160870        2310007B03Rik   0
chr1    107511454       107525597       -2      0       chr1    107500891       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511422       107525600       Serpinb2        0
chr1    107511454       107525597       -2      0       chr1    107511509       107525600       Serpinb2        0
chr1    130882067       130887413       0       -3      chr1    130882073       130887411       Il24    0
chr1    134182466       134189942       2       0       chr1    134182169       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182403       134190031       Chil1   0
chr1    134182466       134189942       2       0       chr1    134182408       134190031       Chil1   0
chr1    135818954       135832956       -2      0       chr1    135818597       135833341       Lad1    0
chr1    149520498       149521157       0       -3      chr1    149410728       149449931       Gm29398 70568
chr1    166166500       166166500       -2      0       chr1    166130459       166166510       Gpa33   0
chr2    5379059         5379374         0       -2      chr2    5293456         5676046         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293456         5714762         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5383533         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5569754         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5601302         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5293458         5715350         Camk1d  0
chr2    5379059         5379374         0       -2      chr2    5298992         5715357         Camk1d  0
chr2    32379110        32381837        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32381854        32381858        0       -2      chr2    32379100        32381915        1110008P14Rik   0
chr2    32384633        32384633        0       -2      chr2    32384636        32387739        Lcn2    3
chr2    32384636        32387773        0       -3      chr2    32384636        32387739        Lcn2    0
chr2    120539390       120563799       0       -2      chr2    120506829       120563831       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563851       Zfp106  0
chr2    120539390       120563799       0       -2      chr2    120507433       120563853       Zfp106  0
chr2    125946785       125975527       2       0       chr2    125859108       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125859228       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866111       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866218       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866219       125983587       Galk2   0
chr2    125946785       125975527       2       0       chr2    125866222       125984298       Galk2   0
chr2    125946785       125975527       2       0       chr2    125920565       125983587       Galk2   0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127633311       127633684       0       -2      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    127634841       127656210       0       -3      chr2    127633225       127656695       Mal     0
chr2    164562717       164562717       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562717       164562717       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164562719       164562720       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562560       164568510       Wfdc2   0
chr2    164568498       164568499       -2      0       chr2    164562715       164568506       Wfdc2   0
chr2    181188007       181188388       -2      0       chr2    181187221       181188506       Ppdpf   0
chr2    181188007       181188388       -2      0       chr2    181187342       181188504       Ppdpf   0

```

After that it was just a matter of aggregating the gene names for each unique region with a quick R script:

```
library(dplyr);
data.tbl <-
    read.delim("annotated.tsv") %>%
    group_by(Chr, Start, End, fwd, rev) %>%
    summarise(Genes = paste(unique(Gene), collapse=","), Distance=min(Distance));

write.table(data.tbl, "collapsed_annotated.tsv",
            sep="\t", row.names=FALSE, quote=FALSE);

```

Collapsed output example:

```
Chr     Start   End     fwd     rev     Genes   Distance
chr1    24613189        24613641        0       3       Col19a1 25655
chr1    24613971        24614812        0       4       Col19a1 26437
chr1    35788729        35789387        -5      0       1110002O04Rik   80208
chr1    93151586        93160225        0       -2      2310007B03Rik   0
chr1    107511454       107525597       -2      0       Serpinb2        0
chr1    130882067       130887413       0       -3      Il24    0
chr1    134182466       134189942       2       0       Chil1   0
chr1    135818954       135832956       -2      0       Lad1    0
chr1    149520498       149521157       0       -3      Gm29398 70568
chr1    166166500       166166500       -2      0       Gpa33   0
chr2    5379059         5379374         0       -2      Camk1d  0
chr2    32379110        32381837        0       -2      1110008P14Rik   0
chr2    32381854        32381858        0       -2      1110008P14Rik   0
chr2    32384633        32384633        0       -2      Lcn2    3
chr2    32384636        32387773        0       -3      Lcn2    0
chr2    120539390       120563799       0       -2      Zfp106  0
chr2    125946785       125975527       2       0       Galk2   0
chr2    127633311       127633684       0       -2      Mal     0
chr2    127634841       127656210       0       -3      Mal     0
chr2    164562717       164562717       -2      0       Wfdc2   0
chr2    164562719       164562720       -2      0       Wfdc2   0
chr2    164568498       164568499       -2      0       Wfdc2   0
chr2    181188007       181188388       -2      0       Ppdpf   0
chrX    73716152        73716152        0       -2      Bcap31  0
chrX    167207094       167209162       0       -3      Tmsb4x  0
chrY    30770844        30772724        3       0       Gm29582 32299

```
> 1 votes

# Answer
Another option is BEDOPS `closest-features`, which efficiently locates nearest annotations, overlapping or otherwise.

```
$ closest-features --closest --delim '\t' intervals.bed genes.bed > answer.bed

```

It's easy to reformat output by piping to `awk`, say.

If files are very large, you can parallelize work into "map-reduce" form by specifying the chromosome name as an option:

```
$ closest-features --chrom chr1 --closest --delim '\t' intervals.bed genes.bed > answer.chr1.bed
$ closest-features --chrom chr2 --closest --delim '\t' intervals.bed genes.bed > answer.chr2.bed
...

```

To reduce, just `bedops -u` all the `answer.*.bed` files.

Running these per-chromosome jobs in parallel makes the work take about as long as searching `chr1`, which tends to be the largest chromosome.
> 1 votes

---
Tags: genome, sequence-annotation, bed
---


# Question
Title: STAR-long parameters for aligning RNA ONT reads to genome
Are there any suggested parameters to align ONT reads to the reference genome using STAR-long? For now, I used the parameters suggested here, but I noticed a weird behaviour.

I have RNA reads (*D. melanogaster*) from `R7` and `R9` flowcells, separately. I only selected to analyze `2D` reads in `pass` category.

I have, respectively, 113249 reads for `R7` and 40318 reads for `R9`. I aligned those reads and get (only!) 150 uniquely mapped reads for `R7` data and 8017 uniquely mapped reads for `R9` data. I tried to run again the same command on a different server with fresh compilation, but the output file is consistent with these 150 reads.

However, if I align the same with GMAP, I get 78016 uniquely mapped reads for `R7` and 33523 uniquely mapped reads for `R9`, so I suspect that something went wrong in the alignment run.

I am aware that the two mappers behave very differently, STAR-long being more precise and preferring to report mappings of fewer reads but at better loci, and GMAP being overall less precise, trying to map the most of the reads but at not-so-good loci.

I was wondering if some of you had experience with this and could suggest me the best parameters for RNA reads from ONT?

# Answer
I've had great results using minimap2, particularly when combined with a pre-treatment of Canu for error correction (using *minimap2* for the read-to-read mapping):

```
# correct reads
~/install/canu/canu-1.6/Linux-amd64/bin/canu overlapper=minimap \
   genomeSize=100M minReadLength=100 minOverlapLength=30 -correct \
   -p 4T1_BC06 -d 4T1_BC06 \
   -nanopore-raw workspace/pass/barcode06/fastq_runid_*.fastq
# align reads
~/install/minimap2/minimap2 -p 10 \
   -a ~/db/fasta/mmus/ucsc/mmus_ucsc_all_cdna.idx \
   -x splice <(pv all/4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  samtools sort > 4T1_BC06_all_vs_mmusAll.bam

```

Update: the most recent nanopore basecaller, combined with the most recent minimap2, seems to now do a decent job with mapping, and pre-correction of reads no longer seems necessary. More recently I've been using LAST to map to the transcriptome, and minimap2 to map to the genome. Minimap2 has the ability to use a homopolymer-compressed genome index, which means that the most common consistent error for nanopore reads (i.e. misjudgement of homopolymer length) will not influence the mapping rate.
> 5 votes

---
Tags: nanopore, rna-alignment, star, gmap
---


title: Gene expression data in GC-RMA format?

text:
# Question
Title: Gene expression data in GC-RMA format?
I'm having a hard time retrieving the most expressed genes in a cell line. One of the databases I've found has values for a lot of genes in 60 different cell lines and it seems the data is obtained using GC-RMA. The data is available here in the file *NCI60\_U133A\_20070815.raw.csv*.

My problem is what does the values actually represent? Would it be correct to assume the genes with high values are the ones most expressed in the corresponding cell line?

I've found a short explanation on why there is no unit definition on the page here.

Here is a snapshot of the file:

```
-            huh-7    huh-7      HEK293  ...
1007_s_at    31.727    55.551    360.731
1053_at     402.633   441.622    557.782
117_at        9.968    10.178      9.565
...

```

# Answer
These values come from a microarray so they represent the number of sequences that where attached to the microarray's probes (1007\_s\_at, ...).

In general yes, a higher value means more expression. The normalizations (GC-RMA) applied are to correct for several effects. They help for testing the probability of the given probe to be a false positive.
> 2 votes

---
Tags: database, gene
---


title: Mapping Information from SAM/BAM file

text:
# Question
Title: Mapping Information from SAM/BAM file
I mapped raw illumina reads to longer pacbio reads and I would like to know the following information from my mapping file (SAM/BAM)

* How many PacBio reads are mapped to at least one illumina reads
* A Table showing each pacbio read, the number of illumina reads that mapped, and total coverage (i.e. percentage of mapped regions)

# Answer
> 1 votes
If you sort and index your BAM file:

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

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

```

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

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

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

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

```

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

# Answer
> 0 votes
This infomation can be obtained with

```
samtools index mybam.bam
samtools idxstats mybam.bam

```

You'll get a table with one row per PacBio read, the length, the number of mapped reads aligned to it and the number of unmapped reads aligned to it.

The number of PacBio reads mapped to at least one illumina read is the number of rows where the 3rd column is greater than 0

```
samtools idxstats mybam.bam | awk '$3>0' | wc -l

```

---
Tags: read-mapping, samtools, sam
---


# Question
Title: Removing repeated reads from nanopore 1D² reads
*This is summary of discussions and a question that was posted by Devon O'Rourke on Twitter*

Following Albacore basecalls on a 1D² library I get two sets of .fq files, summary stats, etc.: one for the 1D basecalling script, and one for the 1D² script. I'd like to use the reads generated from these scripts in a genome assembly.

After a little bit of `grep` searching, it seems like there are overlapping reads among the 1D and 1D² directories - namely, those 1D reads that generated 1D² reads!

How do people approach the choice of selecting which directory to pull the unique read from if you're planning on assembling a genome? I'd think that the 1D² reads are of higher quality and be preferred, but perhaps there are 1D² reads which are lower in quality than one of the two 1D pairs?

Just looking for advice on how folks approach 1D² libraries, fastq outputs, and subsequent mapping. For what it's worth I'm hoping to use Minimap2 to proceed with alignments but am happy to hear others thoughts on alternatives.

## Update

I tried the first solution proposed by @gringer, but my filtered\_1d\_reads.fq file contains same number of reads as raw\_1d\_reads.fq file.

I'm pretty sure the problem is the naming convention difference between the linear and paired reads. linear reads always follow a read\_id of alphanumeric (8)-(4)-(4)-(4)-(12), while paired reads are (8)-(4)-(4)-(4)-(20)-(4)-(4)-(4)-(12). When I looked at it in more detail, the paired read "read\_id" naming is just a concatenation of the two linear "read\_id"'s that are partnered. So that (20) value is really the last (12) of the first linear read, and the first (8) of the second linear pair.

# Answer
You can exclude the 1D² reads by generate a list of the ids from the 1D² reads, then filtering those reads out from the 1D fastq files. Here's one way to do that using some accessory scripts that I have written:

```
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | awk '{print $2}' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq

```

In order to deal with the concatenated read IDs, it's necessary to switch to a language like Perl that supports the use of backreferences for the splitting. The following perl command (replacing the first line in the above procedure) inserts a line break after the 36th character in the 2nd value, then prints that value:

```
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt

```

So a full sequence of operations would look something like this:

```
$ cat 1d2_reads.fq | ~/scripts/fastx-length.pl | perl -lane '$F[1] =~ s/(.{36})/$1\n/; print $F[1];' > readNames.txt
$ cat 1d_reads.fq | ~/scripts/fastx-fetch.pl -v -i readNames.txt > filtered_1d_reads.fq
$ cat 1d2_reads.fq filtered_1d_reads.fq > unique_reads.fq

```
> 1 votes

---
Tags: nanopore, 1d2-reads
---


title: Looper: Duplicate sample name error

text:
# Question
Title: Looper: Duplicate sample name error
I am using Looper. While running `looper run project_config.yaml` I am getting `duplicate sample name` error. I have two `fastq` files for each sample, so in the `annotation.csv` file I have:

```
sample_name,library,file
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r1.fq.gz
abeta-24-1,RNA-seq,data/merged/abeta_24h_1_r2.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r1.fq.gz
abeta-24-2,RNA-seq,data/merged/abeta_24h_2_r2.fq.gz

```

My `pipeline_interface.yaml` file looks the following way:

```
protocol_mapping:
  RNA-seq: sandro_rna_seq

pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-L": "IU"
      "-T": "gencode_mouse_m13"
      "-E": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/jyvo_experiment-metadata.yaml"
      "-C": "/scratch/nv4e/binf-
pipelines/bulk_rna_seq/salmon/pypiper/salmon_rna_seq.yaml"
      "-X": "True"
      "-P": "6"
      "-I": sample_name
    resources:
      default:
        file_size: "0"
        cores: "20"
        mem: "120000"

```

Why am I getting this error? How could I fix it?

# Answer
As explained in the pep docs, the `sample_name` column should be considered a unique identifier for a sample. The PEP requires that **One row corresponds to one sample** You have multiple rows with the same name; this is not allowed.

You simply need to recast your annotation sheet such that it has `file1` and `file2` columns, but only 2 rows.
> 2 votes

---
Tags: scrnaseq, looper, yaml
---


# Question
Title: predict the foldability of single-stranded DNA molecules
I have a list of regions of the human genome and I want to predict if single-stranded molecules in a buffer would tend to fold and create pin structures by sequence self-complementarity. What's the most precise software/parameter set I can use to such effect?

TRIED SO FAR: Vienna RNALFold (unsure if applying an RNA software to DNA is adequate):

```
Compute locally stable RNA secondary structure with a maximal base pair span.
For a sequence of length n and a base pair span of L the algorithm uses only
O(n+L*L) memory and O(n*L*L) CPU time. Thus it is practical to "scan" very
large genomes for short RNA structures.
Output consists of a list of secondary structure components of size <= L, one
entry per line. Each output line contains the predicted local structure its
energy in kcal/mol and the starting position of the local structure.

```

Example:

```
samtools faidx $ref $region | seqtk seq -A - | $HOME/viennarna/bin/RNALfold --noconv
>chr1:100000-100150
.((((.(......).)))). ( -2.70)  127
.(((.((........)).))). ( -5.70)  123
.(((((((.((........)).)))..)))). ( -6.30)  119
.(((.(.(((.((.....))))).)..))). ( -7.00)  117
.((((.((........)).)))). ( -8.30)  109
.((((..((.((..(..(((.....)))..)..)))))))). (-10.10)  103
.((((.((..((..(((...((((.((........)).)))).....)))))..)))))). (-10.50)   90
.((.((.((((....................)))).)).)). ( -4.10)   77
.((((.((.((((....................)))).)).)))). ( -5.70)   75
.((((..((((.((....)).)))).)))). ( -2.40)   72
.((.((..((((..((((.((....)).)))).))))..)))). ( -5.50)   65
.(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).). (-16.20)   40
.((((......(.((..(..................((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).)))).)..)).).....)))). (-18.10)   30
.(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......))). (-18.50)   23
.((....(((....(((((((((...(.........)..)))).)))))((.((..((((..((((.((....)).)))).))))..))))..((((.((........)).))))......)))....)). (-18.70)   17
.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))). ( -9.60)    7
.((.(((((((((......((...(((((..((((...((...))....))))..)))))...))......))))))))).)). (-10.90)    4
.(((..(((((........)))))..))). ( -3.60)    1
ACTTAAGTTGTAGAAGGGAATAACGCAAGAGTGAATTTAGGGCGGGGCAAAAGGATAAATTTTACGGTACAAAGTTTCTACGGGTTTATATATGTATAACTAAGTCCAAGCGCGGGGGATATGGCCAGTGCACAACGGCGGGCATCATAAT
 (-21.70)

```

# Answer
> 1 votes
You could also take a look at Infernal.

I think your best bet is to use RNA-focused programs; I'm not aware of a DNA-specific one, which makes some sense because DNA is typically double-stranded, so the secondary structure community is primarily focused on RNA. I would predict that secondary structure potential of single-stranded DNA would be similar to that of RNA.

---
Tags: phylogenetics, genomics, structural-biology
---


title: Looper: missing argument error

text:
# Question
Title: Looper: missing argument error
While trying to run looper I am getting the following error:

> Error (missing attribute): 'sandro\_rna\_seq' requires sample attribute '/scratch/nv4e/sandro\_looper/pipeline/jyvo\_experiment-metadata.yaml' for option/argument '-E'
> Not submitted: Required attribute(s) missing for pipeline arguments string

In my `pipeline_interface.yaml` file I have the following:

```
arguments:
  "-L": "IU"
  "-T": "gencode_mouse_m13"
  "-E": "/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-
metadata.yaml"

```

The file `"/scratch/nv4e/sandro_looper/pipeline/jyvo_experiment-metadata.yaml"` is definitely present but it still gives an error...

Any suggestions would be greatly appreciated.

# Answer
You are using the argument values incorrectly. Argument values are sample attributes, not full character paths. Hence, "IU" is a sample attribute, "gencode\_mouse\_m13" is a sample attribute, so why would you think "/scratch..." would be treated any differently?

The pipeline interface maps sample attributes to command-line arguments. It is modulated by sample (different values per sample). What you're trying to do is pass a constant string to an argument.

For this, there is a `pipeline_args` section of your project config. Please read the docs on project config.
> 2 votes

---
Tags: scrnaseq, yaml, looper
---


# Question
Title: Splitting characters of a column in repeated lines and reducing them to one line
I have a file with sample data like this:

```
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042017    2   5008    T:5006  C:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     ß@  5
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042017    2   5008    T:5006  C:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     ßß  2
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042030    2   5008    C:5007  A:1     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     @@  1
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042030    2   5008    C:5007  A:1     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     @ß  3
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042029    2   5008    G:5006  T:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184337    2   5008    C:5007  T:1     @@  1
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90042029    2   5008    G:5006  T:2     90184334    90184354    GAGGCCACAAGAGGGCACAA    90184336    2   5008    A:5007  G:1     @ß  3
1   94853396    94853416    GCCGCCACTAGATGGTGCTA    94853398    2   5008    T:5007  C:1     94969254    94969274    CAGTCCTCTAGAGGGAGCCC    94969266    2   5008    G:5006  A:2     ß@  5
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584300   2   5008    C:5006  T:2     103841436   103841456   CATGCCACAAGAGGGCATCA    103841456   2   5008    A:5006  G:2     @ß  3
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584300   2   5008    C:5006  T:2     103841436   103841456   CATGCCACAAGAGGGCATCA    103841445   2   5008    C:5006  T:2     @@  1
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584284   2   5008    A:5007  G:1     103841436   103841456   CATGCCACAAGAGGGCATCA    103841456   2   5008    A:5006  G:2     ßß  2
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103584284   2   5008    A:5007  G:1     103841436   103841456   CATGCCACAAGAGGGCATCA    103841445   2   5008    C:5006  T:2     ß@  5

```

As seen, some lines have repeating values for `$2` and `$3`, line 1-6 and line 8-11 in this case.

`$20` of my file contains codes of two characters representing each line, in the format "@@","ßß", "@ß","@N","ß@","ßN","N@","Nß" or "NN". Each line can have one of these 9 possible codes.

I want to write a script, which splits the characters of `$20` (@,ß or N) of the lines which are repeated (1-6 & 8-11 here), and reduces them into one line in the following manner:

Count the number of @, ß or N in the first and second positions respectively of the repeated lines. For the first and second characters of each repeated line, if:

1. @ occurs most, print @
2. ß occurs most, print ß
3. N occurs most, print N
4. @ occurs equal to ß, print N
5. @ occurs equal to N, print @
6. ß occurs equal to N, print ß

In this sample dataset with repeated lines, the output file generated should be something like this, only columns `$1`, `$2`, `$3`, `$4`, `$11`,`$12`,`$13`, and a new column with the new "@","ß","N" codes after said counting is reported:

```
1   90042011    90042031    AGTGCCACCAGGTGGGGCCG    90184334    90184354    GAGGCCACAAGAGGGCACAA    @N
1   94853396    94853416    GCCGCCACTAGATGGTGCTA    94969254    94969274    CAGTCCTCTAGAGGGAGCCC    ß@
2   103584283   103584303   TGGGCCACTAGGAGGCACTG    103841436   103841456   CATGCCACAAGAGGGCATCA    NN

```

In the first repeated case(1-6), for the first character @ occurs four times and ß occurs twice and for the second character @ and ß occur equal number of times, so the final code is @N. In the second repeated case(8-11), @ and ß for both character positions occur equally, so the final code is NN.

# Answer
> 2 votes
OK. This is a bit convoluted, but it works. There might be simpler solutions...

First we'll write a function that takes a data frame and a character position (ie; first or second position in the string), groups the data frame by the second and third columns, counts occurrences of each character in the specified position of the 18th column, and returns a vector according to your ranking scheme (this doesn't include counts of N since it doesn't occur in your example data, but it's a simple extension of the `case_when` line to add it):

```
library(dplyr)
library(tidyr)
library(stringr)

summarise_characters <- function(df, pos) {
  df %>% mutate(character=str_sub(X18, pos, pos)) %>%
  group_by(X2, X3, character) %>%
  summarise(n=n()) %>%
  ungroup() %>%
  spread(character, n, fill=0) %>%
  mutate(pos1=case_when(`@` > ß ~ '@',
                        `@` == ß ~ 'N',
                        `@` < ß ~ 'ß')) %>%
  `$`(pos1)
}

```

Next we remove duplicate lines from the data frame, then create a new column, pasting together the output of our function called on the first and second character:

```
df %>% group_by(X2, X3) %>%
  slice(1) %>%
  ungroup() %>%
  mutate(counts=paste0(summarise_characters(df, 1), summarise_characters(df, 2))) %>%
  select(X1, X2, X3, X4, X11, X12, X13, counts)

# A tibble: 3 x 8
     X1        X2        X3                   X4       X11                  X12       X13 counts
  <int>     <int>     <int>                <chr>     <int>                <chr>     <int>  <chr>
1     1  90042011  90042031 AGTGCCACCAGGTGGGGCCG  90184354 GAGGCCACAAGAGGGCACAA  90184337     @N
2     1  94853396  94853416 GCCGCCACTAGATGGTGCTA  94969274 CAGTCCTCTAGAGGGAGCCC  94969266     ß@
3     2 103584283 103584303 TGGGCCACTAGGAGGCACTG 103841456 CATGCCACAAGAGGGCATCA 103841456     NN

```

---
Tags: r, bioawk, text-processing
---


# Question
Title: Converting Ensembl Gene IDs to Entrez Gene IDs through biomart
Well, I'm trying to convert a list of Human Gene referenced by Ensembl Gene IDs to Entrez Gene IDs. I have been advised to use biomart.

I tried to get a kind of conversion table for all human genes. I don't know if my settings are wrong, but I didn't find any checkbox for Entrez Gene Id in the section `Attributes` \> `External References`. I've just found the Entrez transcript name ID (see below), but it's not what I need.

How can I use *biomart* for such a conversion?

**Edit:**

According to this video tutorial EntrezGene ID should be an option in biomart (see screenshot below).

That video tutorial is outdated.

# Answer
> 8 votes
The ID you need is the NCBI gene ID, which is the same as the EntrezGene ID.

# Answer
> 8 votes
Conversion using `R`:

```
library(biomaRt)

mart <- useDataset("hsapiens_gene_ensembl", useMart("ensembl"))

    genes <- getBM(
      filters="ensembl_gene_id",
      attributes=c("ensembl_gene_id", "entrezgene_id"),
      values=ensembl.genes,
      mart=mart)

```

Where `ensembl.genes` is a vector of Ensembl gene IDs.

---
Tags: gene, public-databases, biomart, conversion
---


# Question
Title: How to retrieve SNPs data of different humans?
How can I retrieve only SNPs from the dbSNP's FTP site. The files are gzipped and organized by human chromosome (i.e, chr\[1-22,X,Y\]), are these for one human? I need only SNPs data for various humans.

# Answer
If you mean that you want snps from individuals, instead of all together, you can find them in 1000 genomes. Here different individuals from different populations are sequenced and variants are called, including snps.

Maybe you should add to your question what your goal is, what would you like to find?
> 4 votes

---
Tags: genome, snp, phylogenetics, human
---


# Question
Title: Looper: can't specify the individual path for one pipeline run
I have multiple samples that I need to run one pipeline on. I am using Looper: https://github.com/pepkit/looper

I have a real trouble specifying a separate location for each pipeline output files. The first time when I did not care about it thinking that `Looper` would just automatically do that, it lead to multiple files being overwritten. Then, trying to avoid it I tried (and it did not work):

> 1. Adding `{sample_name}` to the `output_dir` of the `metadata` section of `project_config.yaml` file:

```
 metadata:
  sample_annotation: sandro_annotation.csv
  output_dir: /scratch/nv4e/sandro_looper_results/{sample_name}
  pipelines_dir: [../pipeline/pipeline_interface.yaml]

```

> 2. Creating a `derived_column` in the `.csv` annotation file:

```
sample_name,library,sample_name_1_read,sample_name_2_read,output_folder
abeta_24h_1,RNA-seq,read1,read2,out
abeta_24h_2,RNA-seq,read1,read2,out
abeta_24h_3,RNA-seq,read1,read2,out

```

`pipeline_interface.yaml`:

```
...
pipelines:
  sandro_rna_seq:
    name: sandro_rna_seq
    path: sandro_rna_seq.py  # relative to this
pipeline_interface.yaml file
    looper_args: True
    arguments:
      "-S": sample_name
      "--input": sample_name_1_read
      "-O": output_folder
...

```

and `project.config.yaml`:

```
...
derived_columns: [sample_name_1_read,sample_name_2_read,output_folder]
data_sources:
  read1: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r1.fq.gz
  read2: /scratch/nv4e/sandro_looper/data/merged/{sample_name}_r2.fq.gz
  out: /scratch/nv4e/sandro_looper_results/results_pipeline/{sample_name}
...

```

In this particular case the parameter `-O` was formed correctly which I was able to see in `sandro_looper_results/submissions/` `.sub` file, however, it was supplied the second time, apparently from the `output_dir` section of the `project_config.yaml` and was overriding the first `-O` command. So, obviously I tried to delete `output_dir` line but it is a required setting for `Looper`, so it can not run without it.

> 3. I have also tried supplying `-O` command in `project_config.yaml` in the `pipeline_args:` section with no success.

I really do not know what I could do at this point. Any suggestions would be greatly appreciated,

# Answer
> 1 votes
Looper does not generate any of these per-sample outputs.

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

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

For real-life example, see this pipeline

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

---
Tags: scrnaseq, yaml, looper
---


title: knnImputation for missing data

text:
# Question
Title: knnImputation for missing data
I have a dataframe containing missing values. I want to impute them with knnImputation as shown below.

Now my question is: are data imputed row wise or column wise and is there a way to give this option in knnImputation?

```
A <-c(1,2,3,4,5, NA, NA)
B <-c(2,3,NA, 2,4, 7, 6)
C <-c(2,3,4, NA, 5, 1,7)
D<-21:27
data<-data.frame(A,B, C, D)
new_data<-knnImputation(data, k = 3, scale = T, meth = "weighAvg", distData = NULL)

```

# Answer
> 4 votes
The data are imputed using values from the same column

Here's a simplified example to show what's going on:

```
A <-c(1,2,3,4,5)
B <-c(2,3,NA, 2,4)
C <- c(2,3,4, 6, 1)
data<-data.frame(A,B,C)
data

#>  A  B C
#>1 1  2 2
#>2 2  3 3
#>3 3 NA 4
#>4 4  2 6
#>5 5  4 1

# I'm using scale = F to make things simpler
knnImputation(data, k = 3, scale = F, meth = "weighAvg", distData = NULL)

#>  A        B C
#>1 1 2.000000 2
#>2 2 3.000000 3
#>3 3 2.594272 4
#>4 4 2.000000 6
#>5 5 4.000000 1

# calculate euclidean distances, dropping column of interest
dist(data[,-2])

#>1        2        3        4
#>2 1.414214
#>3 2.828427 1.414214
#>4 5.000000 3.605551 2.236068
#>5 4.123106 3.605551 3.605551 5.099020

# use k rows with smallest distance to row of interest to calculate weighted average
# using exp(-dist) as weights
weight_sum <- weight_sum = exp(-2.828427) + exp(-1.414214) + exp(-2.236068)
2*exp(-2.828427)/weight_sum + 3*exp(-1.414214)/weight_sum + 2*exp(-2.236068)/weight_sum

#>2.594271

```

If you want to impute using values from the same row, all you have to do is transpose the matrix:

```
t(knnImputation(t(data), k = 2, scale = F, meth = "weighAvg", distData = NULL))

#>     A        B C
#>[1,] 1 2.000000 2
#>[2,] 2 3.000000 3
#>[3,] 3 3.086729 4
#>[4,] 4 2.000000 6
#>[5,] 5 4.000000 1

```

(I had to lower k because there are only 2 extra columns)

---
Tags: r
---


# Question
Title: Mapping bisulfite reads to a reduced size genome
Is it worth mapping bisulfite reads (WGBS) to a reduced size genome?

I'm interested only in modifications in CpG context, thus instead of mapping to a whole genome (human genome) I would:

1. Extract DNA sequence around each CpG (-1000/+1000 bp)
2. Map reads (using bismark) to this subset

My idea is to reduce mapping time, however afraid that there might be some mapping bias or other issues.

# Answer
> 3 votes
CG sequences are pretty common; I expect you'll hit almost all of the genome when looking at a 2kb window centred on every CG, so there wouldn't be much point in doing that for efficiency purposes.

dariober has pointed out that reads coming from elsewhere in the genome could be incorrectly mapped to the reduced genome with high mapq score since the best, genuine genomic region is not present. So even if you are interested in few small regions, you should still map to the entire genome.

# Answer
> 4 votes
Don't restrict the portion of the genome you're mapping to, you'll just increase false-positive alignments by doing so. If you're worried about mapping time, then (A) use more cores or (B) use a different mapper (BSMAP was the quickest out there last I looked, though in my opinion it's fine to wait a bit longer for better quality results (I normally suggest bwameth to people)).

As pointed out by user 719016 in the comments, the above does not apply to amplicon sequencing. With that, you'll get VASTLY better performance if align to a reduced-size genome...and without an appreciably increased false-positive rate (unless your target-enrichment for sequencing is poor).

---
Tags: ngs, read-mapping, sequence-alignment
---


# Question
Title: remaining human genome variation that hasn't been sequenced?
Given all the genome variation information we have accumulated from resequencing human genomes up until now (early 2018), how much variation is left that hasn't been sequenced?

EDIT: given one of the answers, and the fact that 99% of the SNPs at 1% or more are found, is there a quantitative measure of remaining variation?

# Answer
> 5 votes
In theory, almost any base in the human genome may mutate, so you have billions of variants to go. Ok, this is not so useful.

A related and potentially useful question is: given a human, what is the fraction of his/her variants seen in the the previously sequenced samples? For a Wright-Fisher population, there is an analytical answer: x% of variants on a haplotype has frequency below x%. Replace "x" with "1": 1% of variants on a haplotype have frequency below 1%. The 1000 genomes project concluded that 99% of SNPs at frequency 1% or higher have been called by the project. If we all came from a Wright-Fisher population, 99% of of SNPs have been called in 1000g. Given that the project also calls SNPs at lower frequency, nearly all SNPs in a newly sequenced sample are present in 1000g.

The above is a theoretical analysis. In practice, there are two major complications. First, we are not a Wright-Fisher population. Due to recent population expansion, each haplotype harbors more new mutations than the Wright-Fisher model would predict. Given a known population history, it is actually possible to numerically compute $f(x)$, the fraction of variants with frequency below x (for Wright-Fisher, $f(x)=x$). I don't have this result, though. Second, 1000g and many other genome projects didn't call SNPs in repetitive regions and missed a significant portion of high-frequency indels and structural variations. With Illumina, you can only call ~30% of long deletions callable with PacBio data, but we only have a dozen of PacBio human genomes in public. There are still a lot we don't know.

---
Tags: genome, snp, human, variation
---


# Question
Title: Intron retention events for a condition sample with 3 replicates
I have 3 biological replicates of RNA seq data for a particular condition. I want to find out intron retention events from those biological replicates for a given condition. There is no comparison I am doing here. Just from the three replicates, I am trying to find out what introns are retained in my sample. There is only one condition here. Tools like rMATS require 2 conditions. Is there any tool I can use to find out intron retention events given that I only have one condition (with replicates).

Ideally I would like to use a method similar to that used in this paper, though I can't fully grasp their methodology. This would allow me to get p-values for enrichment of specific introns.

# Answer
> 1 votes
All introns will have some level of retention and I guess no intron will be 100% retained. In the absence of a comparator, you will need to work out how much intron retention an intron needs to be called retained. This is a human judgement, there is no automated way of deciding this.

What you need to do this is some sort of PSI calculation (percent-spliced-in). MISO will provide this, but it doesn't do replicates, so I guess you'll need to run it on each sample and then average the PSIs and determine if it surpasses some sort of threshold you have set.

# Answer
> 0 votes
The method you linked to is interesting, though I suspect that Ian's answer is both easier and generally sufficient to get what you want.

Having said that, the general idea behind the method you linked to is to use DESeq2 to test for differentially expressed introns. To generate the control samples, they're effectively taking the number of reads across all introns in a gene and distributing them over the introns. They weight the distribution such that longer introns receive more reads (they use an effective length for this, but I'm skeptical that this makes that much of a difference). DESeq2 is then giving the p-values you're after. Unless you can get the authors to send you there code, you'll have to write a bit of code to do all of this...but it shouldn't be too bad.

BTW, you can use the `GenomicFeatures` package from bioconductor, specifically the `exonsBy(..., by="gene)` follows by `gaps()`.

# Answer
> 0 votes
IRFinder can be used to detect the IR events and can easily reach your goal.

Here is an example. IRFinder can calculate the IR ratio to detect IR event. mimirna.centenary.org.au/irfinder/example1.html. (You just need run IRFinder for one condition and each replicate will output a folder). Two conditions can be used to detect the differential IR events. If you have only one condition, you can select the IR events from the common events of the three replicates or select IR events exist at least two of three replicates based on the main output files( named IRFinder-IR-dir.txt). The IR events can be selected by filtering the IR ratio, Intron coverage etc

---
Tags: rna-seq, rna-splicing
---


# Question
Title: Write the output of the VariantAnnotation package to file
I'm having some problems in exporting the annotations of some variants from a VCF file, after using the "VariantAnnotation" package from bioconductor. Below the script I've used:

```
#!/usr/bin/env Rscript

library(VariantAnnotation)
library(BSgenome.Scerevisiae.UCSC.sacCer3)
library(TxDb.Scerevisiae.UCSC.sacCer3.sgdGene)

args = commandArgs(trailingOnly=TRUE)

vcf <- readVcf(args[1], "sacCer3")

sc <- BSgenome.Scerevisiae.UCSC.sacCer3
txdb <- TxDb.Scerevisiae.UCSC.sacCer3.sgdGene

coding <- data.frame(predictCoding(vcf, txdb, seqSource=sc))

```

The `coding` dataframe looks like this:

```
  seqnames start end width strand paramRangeID REF     ALT QUAL FILTER
1     chrI   455 455     1      +         <NA>   C       T   40   PASS
2     chrI   455 455     1      +         <NA>   C           40   PASS
3     chrI   455 455     1      +         <NA>   C       G   40   PASS
4     chrI   455 455     1      +         <NA>   C       A   40   PASS
5     chrI   456 456     1      +         <NA>   T       A   40   PASS
6     chrI   456 456     1      +         <NA>   C CCTGCTC   40   PASS
  varAllele CDSLOC.start CDSLOC.end CDSLOC.width PROTEINLOC QUERYID TXID CDSID
1         T          121        121            1         41       1    1     1
2                    121        121            1         41       2    1     1
3         G          121        121            1         41       3    1     1
4         A          121        121            1         41       4    1     1
5         A          122        122            1         41       5    1     1
6   CCTGCTC          122        122            1         41       6    1     1
   GENEID    CONSEQUENCE REFCODON  VARCODON REFAA VARAA
1 YAL069W  nonsynonymous      CTC       TTC     L     F
2 YAL069W not translated      CTC        TC
3 YAL069W  nonsynonymous      CTC       GTC     L     V
4 YAL069W  nonsynonymous      CTC       ATC     L     I
5 YAL069W  nonsynonymous      CTC       CAC     L     H
6 YAL069W  nonsynonymous      CTC CCCTGCTCC     L   PCS

```

But if I try to save it to file, for instance using the following command:

```
write.table(coding, 'test.tsv', sep='\t', row.names=F)

```

The output table looks like this:

```
"seqnames"      "start" "end"   "width" "strand"        "paramRangeID"  "REF"   "ALT"   "QUAL"  "FILTER"        "varAllele"     "CDSLOC.start"  "CDSLOC.end"    "CDSLOC.width"  "PROTEINLOC"    "QUERYID"       "TXID"  "CDSID" "GENEID"        "CONSEQUENCE"   "REFCODON"      "VARCODON"      "REFAA" "VARAA"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  "T"     121     121     1       41      1       "1"     1       "YAL069W"       "nonsynonymous" "CTC"   "TTC"   "L"     "F"
"chrI"  455     455     1       "+"     NA      "C"     <S4 object of class "DNAStringSet">     40      "PASS"  ""      121     121     1       41      2       "1"     1       "YAL069W"       "not translated" "CTC"   "TC"    ""      ""

```

Notice the `<S4 object of class "DNAStringSet">` bit. Is there any way to force R to give me the same representation as `head`? I couldn't find anything in Biostrings' documentation.

# Answer
> Is there any way to force R to give me the same representation as `head`?

`head` just cuts off everything after the first six entries. The magic, as it were, happens in the default `print` method for the specific type. But the representation generated by it is only useful for ad-hoc display, not for storage in a structured data file.

The best solution is therefore to manually convert the column types. Since it’s only one column:

```
coding$ALT = sapply(coding$ALT, as.character)

```

I’m a bit puzzled that the `REF` and `ALT` column have different types, anyway: that’s probably a bug in the {VariantAnnotation} package that should be reported.

Furthermore, you probably want to use `quote = FALSE` for the `write.table`. In general, `write.table` has terrible defaults. `readr::write_tsv` is almost always a better choice.
> 2 votes

---
Tags: r, bioconductor, variation
---


title: Perplexity is too large

text:
# Question
Title: Perplexity is too large
I am trying to run {Rtsne}. On `pca` dataset prepared by using `dudi.pca` from {ade4}.

For pca calculation I selected 40 components while specifying parameters for `dudi.pca`:

```
pca = dudi.pca(counts(sce_normalized), scannf=FALSE, nf=40)

```

{Rtsne} gives me an error: “perplexity is too large”, how could it be fixed?

The code I am trying to run is below:

```
name_of_file <- "dudi_pca.rds"
path_to_file <- paste0("/scratch/nv4e/scRNA-
seq/normalization/new_datasets/",
                   name_of_file, sep="")
dudi_pca <- readRDS(path_to_file)

# run Rtsne with default parameters
rtsne_out <- Rtsne(as.matrix(dudi_pca), pca=FALSE)

```

# Answer
> 2 votes
I was using `dudi_pca` incorrectly. The supplied parameter to `as.matrix()` should have been `dudi_pca$li`:

```
rtsne_out <- Rtsne(as.matrix(dudi_pca$li), pca=FALSE)

```

---
Tags: r, pca
---


# Question
Title: Picard CollectGcBiasMetrics ignoring certain chromosomes/sequences
What's the easiest way to run Picard GCBias ignoring certain chromosomes/sequences in the reference? Looking at the CollectGcBiasMetrics, it seems there isn't a bed file option that can be passed.

```
java -jar ~/picard-tools-1.134/picard.jar CollectGcBiasMetrics --help
USAGE: CollectGcBiasMetrics [options]

Documentation: http://broadinstitute.github.io/picard/command-line-overview.html#CollectGcBiasMetrics

Tool to collect information about GC bias in the reads in a given BAM file. Computes the number of windows (of size specified by WINDOW_SIZE) in the genome at each GC% and counts the number of read starts in each GC bin.  What is output and plotted is the "normalized coverage" in each bin - i.e. the number of reads per window normalized to the average number of reads per window across the whole genome..

Version: 1.134(a7a08c474e4d99346eec7a9956a8fe71943b5d80_1434033355)

Options:

--help
-h                            Displays options specific to this tool.

--stdhelp
-H                            Displays options specific to this tool AND options common to all Picard command line
                              tools.

--version                     Displays program version.

CHART_OUTPUT=File
CHART=File                    The PDF file to render the chart to.  Required.

SUMMARY_OUTPUT=File
S=File                        The text file to write summary metrics to.  Default value: null.

WINDOW_SIZE=Integer           The size of windows on the genome that are used to bin reads.  Default value: 100. This
                              option can be set to 'null' to clear the default value.

MINIMUM_GENOME_FRACTION=DoubleFor summary metrics, exclude GC windows that include less than this fraction of the
                              genome.  Default value: 1.0E-5. This option can be set to 'null' to clear the default
                              value.

IS_BISULFITE_SEQUENCED=Boolean
BS=Boolean                    Whether the SAM or BAM file consists of bisulfite sequenced reads.  Default value: false.
                              This option can be set to 'null' to clear the default value. Possible values: {true,
                              false}

METRIC_ACCUMULATION_LEVEL=MetricAccumulationLevel
LEVEL=MetricAccumulationLevel The level(s) at which to accumulate metrics.  Default value: [ALL_READS]. This option can
                              be set to 'null' to clear the default value. Possible values: {ALL_READS, SAMPLE,
                              LIBRARY, READ_GROUP} This option may be specified 0 or more times. This option can be set
                              to 'null' to clear the default list.

INPUT=File
I=File                        Input SAM or BAM file.  Required.

OUTPUT=File
O=File                        File to write the output to.  Required.

ASSUME_SORTED=Boolean
AS=Boolean                    If true (default), then the sort order in the header file will be ignored.  Default
                              value: true. This option can be set to 'null' to clear the default value. Possible
                              values: {true, false}

STOP_AFTER=Long               Stop after processing N reads, mainly for debugging.  Default value: 0. This option can
                              be set to 'null' to clear the default value.

```

Given this, and wanting to calculate GCBias the using Picard (not other software), what's the best option?

# Answer
The proper solution is to use a different tool (some of this you could do with `computeGCBias` in deepTools). But since you don't want to do that, you'll have to manually remove unwanted chromosomes from the BAM file(s):

```
$ cat foo.awk
BEGIN{split(excludes, excludeList, " ")}
{exclude=0
if($1 == "@SQ") {
    for(ex in excludeList) {
        if(ex == substr($2, 4)) {
            exclude=1
            break
        }
    }
} else { \
    for(ex in excludeList) {
        if(ex == $3) {
            exclude=1
            break
        }
    }
}
if(!exclude) print}

$ excludes="chr1 chr2"
$ samtools view -h alignments.bam | awk -v excludes="$excludes" -f foo.awk | samtools -bo alignments.filtered.bam

```

The awk script is stripping unwanted chromosomes from the header and the body. You can then run the result through picard.
> 1 votes

---
Tags: sam, picard
---


title: TCGA data sets joining multiple files

text:
# Question
Title: TCGA data sets joining multiple files
I have downloaded harmonised the TCGA-HNSC dataset and extracted it. It has 546 samples in total, so what I see in it are FPKM normalised values. The first column remains the same with 60000 rows and 2 columns that is first column ENSEMBL ID and second column FPKM value.

So now the issue is I want to join all the 546 txt file into one where the first column remains the same and rest 546 column would be joined. I'm using this command:

```
paste -d " " *txt > merged_file.txt

```

This is joining the all the data sets as ENSEMBL ID and FPKM value..

But I want 1st column as the ENSEMBLE ID and rest the sample FPKM any suggestion or help would be appreciated

# Answer
> 2 votes
In R:

```
flist = list.files("*.txt", path=".")
fl = lapply(flist, function(x) read.delim(x, header=F))
m = sapply(fl, function(x) x$V2)
row.names(m) = fl[[1]]$V1
colnames(m) = flist

```

That will give you an appropriately formatted matrix named `m` in R, which is probably where you want to continue processing things anyway.

---
Tags: rna-seq
---


title: population admixture

text:
# Question
Title: population admixture
I am new to Admixture and I would like to use fastNGSadmix to look for admixture in one individual (modern). Admixture is supposed to look for proportions of ancestral populations in a sample (either a population or an individual sample). for that I will need to prepare my reference population set to compare with my actual sample which is an individual form region X for instance.

Now I am thinking whether I have to add the X population to my reference set or not? If yes, wouldn't it introduce bias? Like my individual will always show the X population as the first admixing population?

# Answer
> 2 votes
Disclaimer, I am not an expert on admixture, answer is provided without any warranty.

It depends what you want to figure out. Why do you want to categorize the ancestry of the individual? If you analyze this origin of population X you should make a reference without population X and analyze all the individuals from the population with all the surrounding populations in the reference. However if you are interested in origin of a single individual, than the population X got to be in reference as well.

If your result will be that 90% of SNPs of your individual correspond to your population X, then it means that's a normal X fellow. Why it should be a problem?

-- edit --

I just found this preprint that highlights potential pitfalls in interpretation of STRUCTURE/ADMIXTURE results.

---
Tags: reference-genome, phylogenetics
---


# Question
Title: Histogram of different classes of variables in one column
I have a basic ggplot2 R question. I have a dataframe in which data looks like this:

```
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7

```

The file is of multiple lines. I want to make a histogram using gplot2 that will plot the 9 kinds of variables AA,BB,AB,BA,AN,BN,NN,NB,NA on the x-axis and the sum of it's corresponding values in the next column on the y-axis.

So for example, here the y-values should be

```
AA= 13
BB= 12
AB= 5
BA= 4
NN= 7
AN= 10
NA= 3
BN= 12
NB= 1

```

# Answer
This isn't a bioinformatics question, but it's quicker to answer than to close.

```
df1 <- read.table(text = "
AA 4
AA 6
BB 6
AB 5
BA 4
AA 3
NN 2
AN 6
NN 5
AN 4
NA 3
BB 6
BN 5
NB 1
BN 7", na.strings = "")

library(ggplot2)

ggplot(df1, aes(x = V1, y = V2)) +
  geom_bar(stat = "identity")

```
> 2 votes

---
Tags: r, visualization, ggplot2
---


# Question
Title: How to generate a weblogo image that takes into account sequence frequencies?
I have a large aligned amino acid sequences file (fasta file) from next generation sequencing data. The unique sequences have been aligned but these all have certain read counts (frequencies associated with them).

I would like to generate a weblogo image that has the true probability of amino acids at each position. Is there a way to input frequencies/read counts with the aligned files of sequences using weblogo? For example, in the title of the fasta file or as a csv/tsv file?

So that when generating the weblogo image it takes into account the abundances of each sequence?

# Answer
> 3 votes
Don't use the unique sequences. The whole point of sequence logos is that the height of the letter represents its frequency at a given position. You remove this information when you make your sequences unique. So, while it makes sense to use unique sequences when aligning, you shouldn't when building your logo.

Generate a simple mapping file (`counts`) with the number of times each sequence was found:

```
seq1 3
seq2 10
seq3 2

```

Then, given an aligned protein fasta file (`foo.aln`) like this:

```
>seq1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

```

You can simply repeat the sequences as needed (see my answer here for the `FastaToTbl` and `TblToFasta` commands):

```
$ awk '
        (NR==FNR){
            count[$1]=$2;
            next
        }
        {
            for(k=1;k<=count[$1];k++){
                print $1"."k"\t"$2
            }
        }' counts <(FastaToTbl foo.aln) | TblToFasta
>seq1.1
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.2
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq1.3
MEEPQS---DPSVEPP-LSQETFSDLWKLCFLPENNVLSPLPSQAM-DDL
>seq2.1
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.2
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.3
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.4
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.5
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.6
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.7
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.8
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.9
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq2.10
MEEPQSCFGDPSVEPPPLSQETFKDLWKL--LCENNVLS---SQAM-DDL
>seq3.1
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL
>seq3.2
MEEPCS---DPSVEPPQLSQETFSDLRKL--LPENNFLSPLPSQAMCDDL

```

Note that the above assumes your shell supports the `<(command)` syntax for process substitution. Bash, the default shell on most Linux and macOSX installations does, but if you are using a shell that doesn't, you might need to run `FastaToTbl foo.aln > foo.tbl` first, and then use `foo.tbl` as input.

---
Tags: sequence-alignment
---


# Question
Title: What is the standard way to work with a diploid reference genome? Complementary strands?
At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

More and more biologists are using emerging technologies into order to capture the diploid nature of genetic information, e.g. phasing SNPs between mother and father chromosomes.

What is the standard way that bioinformaticians generated a standard diploid reference genome?

Actually, reference genomes are not truly haploid (by my understanding). Given that reference genomes are 5'-3' coordinated, in order to create a complementary strand, one would need to take the 3'-5' complement. In order to have a diploid genome, you need two reference genomes and two 3'-5' complements.

More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

# Answer
For calling small variants, the *standard way* is to simply call diploid genotypes. You can already do a variety of research with unphased genotypes. You may further phase genotypes with imputation, pedigree or with long reads/linked reads, but not many are doing this because phasing is more difficult, may add cost and may not always give you new insight into your data. For these analyses, we use a haploid genome. For human samples, the vast majority of "large-scale genomic studies" are done this way.

A diploid reference actually doesn't help much with reference-based analysis; it only complicates algorithms. What could help a lot is a population reference, which may be represented by a graph or a compressed full-text index or both. In theory, if you have a comprehensive population reference and a capable mapping algorithm, you may call extra variants that would not be callable with short reads. In practice, however, there are quite a few technical challenges. Handling population references is a research topic. There are no "standards" yet.

If the goal is to assemble a new reference genome from a diploid sample, we almost always prefer to produce a diploid assembly. Unfortunately, I believe there are no "standard" procedures, either. SuperNova from 10x genomics builds the diploid information into a graph. Falcon from PacBio uses "unzip". I don't think they have got widely used and evaluated so far.

PS: saw your edit while writing the above. The fact that the genome only represents one strand does not mean we have to create the complement strand explicitly in analyses. We do most of reverse complement on the fly in algorithms as well as in mind.

> reference genomes are not truly haploid

That depends on how the reference is assembled. If you sequence a haploid sample (e.g. bacteria), your assembly will be haploid. If you sequence an inbreed lab strain that is almost homozygous (e.g. mouse and fruit fly), your assembly will be nearly haploid. If you sequence a diploid sample, your assembly is very likely to be a mosaic of the two haplotypes. In case of the human reference genome, it is more complicated. It is largely a mosaic of several humans by stitching ~150kb haplotypes from these samples.
> 5 votes

# Answer
> At the moment, the standard reference genomes (e.g. hg19, hg38) are haploid genomes. We know that the human genome is diploid. Naturally, the latter would be the respectively correct representation of the human genome.

The premise of the OP's question is false. The natural reference representation of the human genome is not diploid.

### Think of a reference genome as a map, and not as a specific example of a human being's DNA.

Not only is the human genome reference haploid, but it is also a composite genome. This means that the human genome reference sequence is composed of sequences from multiple individuals. In other words, the human reference does not correspond to any one human sequence.

Any particular read from a DNA sequencer will be a read from a human genome that diverges from the reference genome. So an algorithm that tries to match the read to the reference genome will always need to handle potential discrepancies. Adding a second map against which to match a read would not change that fact. Therefore, there is little value in providing a second haploid reference genome.

**Side note**: There are parts of the human genome that "are too complex to be represented by a single path", and the Genome Reference Consortium provides "alternate loci" for such regions of the genome.

### Answer to the OP's second question

> More importantly, how have large-scale genomic studies dealt with the fact that the haploid reference genome is a consensus-based "half" of a human genome?

I interpret this question as "How do large genomic studies represent genome diversity that cannot be represented by a haploid reference genome?"

One standard approach that has been very popular in genome-wide association studies is the use of a haplotype reference panel such as for example from the 1000 Genomes Project. Modern phasing programs exploit the shared ancestry of the samples in a study and reference haplotypes to phase a study's samples. One popular phasing program is Impute2
> 6 votes

# Answer
To add to all the other great answers, I would mention that the question is somewhat misleading. If the reference genome is for a single individual, then it should be diploid. However, it's a reference for all humans. It should really contain billions of copies to fully account for all the diversity. Since that is not realistic, the reference serves as a simple approximation.

This point was addressed by the recent Korean genome paper:

> Human genomes are routinely compared against a universal reference. However, this strategy could miss population-specific and personal genomic variations, which may be detected more efficiently using an ethnically relevant or personal reference. ... Systematic comparison of human assemblies shows the importance of assembly quality, suggesting the necessity of new technologies to comprehensively map ethnic and personal genomic structure variations.
> 3 votes

# Answer
There are some assemblers that produce assembly graphs that attempt to describe all the possible haploid paths within a set of reads. Such an assembly attempts to capture all the diploid variation (and/or population variation) in a sample at the expense of not having full-length chromosomes.

Canu (for example) will produce contigs that are extended as long as consensus is maintained across different reads, but when there is a reliable break in coverage (i.e. an area where chromosomes are heterozygous) then the contigs will be broken up. Canu provides as output a GFA file (assembly graph) that can be used to determine which paths might combine together into a single chromosome.
> 2 votes

# Answer
Since nobody has tackled you second question, I'll give it a go:

The complementary strand to a genome sequence is the 5'-3' *reverse* complement, not the 3'-5' complement. This sequence is not captured in standard reference genomes, but the information is. Given a sequence, it's trivial to calculate it's reverse complement, and any tool that's designed to work with sequence data takes this feature into account.
> 1 votes

---
Tags: human-genome, genome, reference-genome
---


# Question
Title: How can I subset a reference based on only the first chromosome?
I have the GATK FASTA reference for hg19. I would like to subset this based on the first chromosome. Normally what I do is open the FASTA as a text file, and then delete all other chromosomes except chromosome 1.

That's awful.

Picard tools has an option. But this doesn't appear to subset only by chromosome. How could this be done efficiently?

# Answer
It sounds like you want `samtools faidx foo.fa` followed by `samtools faidx foo.fa chr1 > your_subset_file.fa` (or whatever the first chromosome is). The output file is then a regular fasta file subset as you like. You can get any chromosome you want with that. In fact, you can also do regions (e.g., `chr1:100-1000`, though note that the sequence name in the resulting fasta output will be `chr1:100-1000`, which may break some downstream tools).
> 3 votes

---
Tags: fasta, reference-genome, picard
---