text:
Writing a perl script to holding information for two genes
Basically I have a perl script in which I have an array (where each element of the array references a hash) and need to be able to print the array with a dumper function.
Thus I need to be able to split the $line on white spaces and save into variables and then construct anonymous hash and push into @genes:
#!/usr/bin/perl
use strict;
use warnings;
use Data::Dumper; # print complex data structure
my @genes; # declare the array
while(<>) { # this means that as long as lines come from the pipe we keep going
  my $line = $_; # a line that come from the pipe (we go line by line)
  next unless $line =~ /^\d+/; # skip lines except those reporting genes
  <split the $line on white spaces and save into variables>
  <anonymous hash and push into @genes>
}
print Dumper(\@genes);
exit;

I don't  even know how to approach this, this is my first time perl scripting and I am really confused. Any insight would be extremely helpful.
So far I have:
my $gene = {'id' => $id, 'start' => $start, 'end' => $end, 'frame' = $frame, 'score' => $score};
my @line_split = split(' ', $line);
my $id = $line_split [0], $start = $line_split [1], $end = $line_split [2], $frame = $line_split [3], $score = $line_split [4]

What is the format of mys.coord2?

00001      16     324  +1   0.929
00002     751     308  -2   0.911

simply just contains two lines

Hi code_pink now the question is more specific, but we still don't know how is it failing, or why. Could you [edit] it further to include this information? Also if you would like to have some guidance explaining your reasoning with your code would be useful to explain why it is a good direction or not (this is a great tool to find the way to code). Also if you mark one answer as accepted, people (usually) won't answer further.

Perhaps this small test script will help demonstrate some of the principles:
#!/usr/bin/env perl

use strict;
use warnings;
use Data::Dumper;

my @arr;

# define the current line
my $line = "foo\tbar\tbaz\n";

# remove the trailing newline character
chomp $line;

# put the line's elements into an anonymous array
my ($a, $b, $c) = split('\t', $line);

# insert key-value pairs into anonymous hash
my $hash = {a => $a, b => $b, c => $c};

# add the hash table to the array
push @arr, $hash;

# print the array
print Dumper \@arr;

Run this script and see what it prints out.

`my ($id, $start, $end, $frame, $score) = split('  ', $line);
my $hash = {id => $id, start => (at) start, end => (at)end, frame => (at) frame, score => (at) score};`

is this sort of in the right direction

my ($id, $start, $end, $frame, $score) = split('\t', $line); my $hash = (id => $id, start => $start, end => $end, frame => $frame, score = $score);

Sorry, that should be: my $hash = { ... } not using parentheses, but brackets. I can't edit my previous comment so hopefully this helps.

The following Perl documentation pages should be informative:

split - for splitting a scalar at all matches of a defined pattern
perlreftut - discusses the approach of anonymous variables and how to combine them


text:
How to efficiently compute the exact percentage of non-unique k-mers in a genome for given k?
I'm looking for some software that can "efficiently" (time and memory) compute the exact percentage of non-unique k-mers in a genome for given k. I don't need the k-mers or the abundances itself, I just need the percentage.
Alternatively, the result could be:

the number of different k-mers observed in the genome
the number of unique k-mers, i.e. k-mers that occur only once in the genome, or the number of non-unique k-mers, i.e., k-mers that occur more than once in the genome.

which can be easily used to compute the requested percentage.
Requested features:

don't count k-mers with ambiguous nucleotides
count a k-mer and its reverse complement together
k upto 50 or 100
for large genomes (upto some GB)
efficiency

For the study in mind, I like to vary k. Hence, I will run the software several times with different k.
Any help is welcome.

Welcome to the site. What is your current solution? Which non-efficient programs have you found (To avoid proposing them)? Can't the other answers in the [tag:k-mer] help you to make that program?

You can do all of that with khmer. For example, abundance-dist-single.py produces a file with columns: k-mer abundance, k-mer count, cumulative count, and fraction of total distinct k-mers. So for question 1 you would sum column 2. For question 2 you would just get thek-mer countassociated with ak-mer abundance` of 1.
That package also provides a python API (see the read the docs link above) if you need to customize things more.

Thanks a lot for that hint. I made a first test of khmer using A.thaliana: python abundance-dist-single.py -k <k> -b <genome> <out>. It looks great. It's fast and easy to use. However, I get an error for k > 32. Is there a way to go further than 32?

I'm not sure it's possible with the command line script, but it appears to be possible in the API.

Thanks for that hint, I'll check that. However, I also found out for some examples (and reading the documentation/papers afterwards ;) that it is not exact. Hence, I currently also use the parameters -N 20 -x 1E7 to avoid too many mistakes. Just in case anyone like to do a similar case study. In addition, I have some problems with some species where khmer states ERROR: abundance distribution is uniformly zero; nothing to report. Please verify that the input files are valid. The files seem to be okay and have been used in other analysis. Any ideas?

No, sorry, I've never seen that error.

source: Stack Exchange

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/

sort: votes


text:
Counting repeated kmers sequences that match at least x % of reads sequence
Working on a fastQ file, I would like to get the occurrences of repeated sequences for all possible kmers of a given length that cover at least 90% of the read's length for the whole data set.
example :

for a length 6 with the kmer "ATTGGG" and a data set  containing reads of 300 bases each

I want to obtain the occurrences of the following kmers :
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGG......  (90 % of the read's length without  gaps)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGA.....  (91%)
ATTGGGATTGGGATTGGGATTGGGATTGGGATTGGGAT....  (92%)

this for all possible combinations of kmers of length = 6.
I have done some research and tested many softwares such as (KMC2,Jellyfish,Khmer,Dsk,Scturtle) but none matched what I would like to do, because I want to obtain the occurrences of the repeated kmer without any possible gaps in-between with a certain % of read coverage.
I've done a naive algorithm that works but the optimization isn't the best obviously. Is there a tool that does the trick ?

The following javascript does what you want. You need node.js to run it. It should be easy to translate the code to Python. I have not carefully tested it. Use with caution.
EDIT (response to new comments): the program has been changed to compute the sum of lengths. Note that only a tandem repeat of length k*2 or longer is counted. For example in sequence ATTGGGATTGGGATTcGGGATTGGG, the script returns 15, because the second is not long enough. Counting length smaller than k*2 will be slower and more complex, but is often not the right thing to do. For example, in ATTGGGATTGGGcccccccgaaatcgatagcatcgaGGGATTGcgatc, we would not count the second GGGATTG as a tandem repeat.
On performance, for an input string of length $l$, the time complexity of this script is $O(l)$ associated with a small constant. It will be efficient and adequate to practical uses.
If you use TRF, you should run it on your sequences and then test if any tandem repeats it finds match GGGAAT – the procedure will be similar to my script, except that TRF will also find impure repeats.
// given the repeat unit length _k_, find all tandem repeats of 2*k or longer
function trf_k(k, str)
{
    var streak = 0, a = [];
    for (var i = k; i <= str.length; ++i) {
        if (i < str.length && str[i] == str[i - k]) {
            ++streak;
        } else {
            if (streak >= k)
                a.push([i - streak - k, i]);
            streak = 0;
        }
    }
    return a;
}
// find the sum of tandem repeat length with _kmer_ being the repeat unit
function trf_kmer(kmer, str)
{
    var a = trf_k(kmer.length, str), sum = 0;
    for (var i = 0; i < a.length; ++i)
        if (str.substring(a[i][0], a[i][1]).indexOf(kmer) >= 0)
            sum += a[i][1] - a[i][0];
    return sum;
}

var str = 'xATTGGGATTGGGATTGxGATTGGGATTGGGATTGGGATTGxx';
console.log(trf_kmer('GGGATT', str));

Sorry i mean that i want the perfect match of each possible repeated sequences, and in the case of your algorithm, when there is a gap it's still count the kmer. Indeed the task is highly customized. I guess I'll just start using parallelism on my program. I'll take a look into TRF, thanks.

@hilta007 Thanks for the explanation (I am deleting my old comment BTW). If you want to find the total length, just change "max" to "sum", but with a caveat (see edits for detail). My implementation should be fairly efficient. I wouldn't worry too much about performance. A good single-threaded script can sometimes be faster than a bad multi-threaded C/C++ program.

Thanks for the edits, waiting the results of TRF and my own algorithm to compare quality/speed for now. I'll keep your algorithm in mind.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/3.0/


text:
Run gffread in multi-thread mode
Is there any possibility to run gffread in multi-thread mode? The answer seems to be 'no' from the manual (or gffread -h), as no multi-thread option is mentioned.
I'm mostly using this utility to extract transcript sequences (FASTA) from annotation files (GTF). In single-thread mode, runtime is acceptable in most cases, but would like to speed up a bit in some cases (especially for big genomes).
Note that I already tried, as mentioned in the manual, to index the genome first (which, again, cuts the runtime in most but not all cases):

Note that the retrieval of the transcript sequences this way is going to be much faster if a fasta index file (genome.fa.fai in this example) is found in the same directory with the genomic fasta file.
Such an index file can be created with the samtools utility prior to running gffread, like this: samtools faidx genome.fa
Then in subsequent runs using the -g option gffread will find this
  FASTA index and use it to speed up the extraction of transcript
  sequences.

I'm also open to different alternatives than gffread to convert GTF -> FASTA.

I have used gffread a few times. My impression is it is fast. As I have just tried, for Ensembl human annotations, it gets all transcripts in 29 sec. This should be faster than many other operations on human genome. How often do you want to run gffread?

Thank you. Yes it is usually quite fast. Actually I was adressing this question since it took me more than 8min to get the transcripts out of a putative genome size of 1.4Gb, but could not reproduce the slow runtime on another server (where it took me around 110sec). Weird.

Since there doesn't seem to be an easy way to run this in parallel, you could instead break the job into sections. For example, separate each chromosome into its own gtf file, extract the sequences using that file and then cat them all together.
The commands would be something like this (using the human hg38 genome and the gff GENCODE annotations downloaded from here):
Split the file into one file per chromosome (or whatever else you have as the first field of your gff/gtf file):
tmp=$(mktemp -d); ## create a temp directory
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf

Extract the sequences described in each file, running each command in the background (&):
for f in $tmp/*; do
    gffread -w "$f.fa" -g hg38.fa "$f" &
done

Concatenate them into one file:
cat "$tmp/*fa" > all.fa

You can combine all this into a single command, using wait to make sure you don't concatenate until all of the sub-commands have finished:
tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' gencode.v27.annotation.gtf &&
    for f in "$tmp"/*; do
        gffread-0.9.12.Linux_x86_64/gffread -w $f.fa -g hg38.fa "$f" &
    done
    wait
    cat "$tmp"/*fa > all.fa

Or, just save it as a script called pgffread.sh (or whatever):
#!/bin/sh
gtf="$1"
genome="$2"

tmp=$(mktemp -d);
awk -vtmp="$tmp" '($1!~/^#/){print > tmp"/"$1".gff"}' "$gtf" &&
    for f in "$tmp"/*; do
      gffread -w "$f".fa -g "$genome" "$f" &
    done
wait
cat "$tmp"/*fa
rm -rf "$tmp"

And then run it giving the genome and annotation files as parameters:
pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa

On the system I ran it on, which has 48 cores, the time difference was significant (although it's so fast anyway, I really don't know if it is worth the effort unless you do this very often):
$ time gffread -w out.fa -g hg38.fa gencode.v27.annotation.gtf
real    1m52.072s
user    1m3.160s
sys     0m47.257s

$ time ./pgffread.sh gencode.v27.annotation.gtf hg38.fa > all.fa
real    0m17.533s
user    1m5.228s
sys     0m43.039s

thanks! yes, I should have thought of breaking the reference into smaller bits.


text:
Using a Bash Script to search TaxIDs against NCBI's Taxonomy yields "400 Bad Request" error?
I've been searching TaxIDs against NCBI's Taxonomy DB to get taxonomic lineages for species.  I have successfully done this for 1,000's of TaxIDs that were returned to me by Blast+ blasts in a CSV.  (The first column in the CSV was the TaxID for the top hit of a sequence, and the other columns were for GI, Scientific Name, e-value, etc.).  I simply wrote a script that would take the first field, which is always the TaxID, and search it against Taxonomy to get lineages for the species those TaxIDs represented.  Everything worked fine!
However, some of the records returned from my original blast (that gave me the TaxIDs in the first place) would return two or more TaxIDs per record.  (I suspect because a sequence matched equally well with two or more different species?)  That's fine, but I had to remove these few exceptions before I searched the rest of the regular results against Taxonomy.  So I've got most of my data/lineages, but for about 1:10,0000 records that gave multiple TaxIDs, I would like to know the lineages of the multiple species that the blast suggested.
I now took all of those "more-than-one-taxid" hits, separated them out into columns in an excel sheet and saved it as a CSV, and am trying to search these against Taxonomy with essentially the same script I used for the rest of the data, but am getting an error for every taxID I search (see below).
Complete Error (for each taxid searched against Taxonomy):
400 Bad Request
&rettype=full&retmode=xml&edirect_os=linux&edirect=7.90&tool=edirect&email=x@xxx.localdomain'my&id=464126
Result of do_post http request is
$VAR1 = bless( {
                 '_protocol' => 'HTTP/1.1',
                 '_content' => '<?xml version="1.0" encoding="UTF-8" ?>
<!DOCTYPE eEfetchResult PUBLIC "-//NLM//DTD efetch 20131226//EN" "https://eutils.ncbi.nlm.nih.gov/eutils/dtd/20131226/efetch.dtd">
<eFetchResult>
 <ERROR>WebEnv parameter is required</ERROR>
</eFetchResult>
',
                 '_rc' => 400,
                 '_headers' => bless( {
                                        'connection' => 'close',
                                        'set-cookie' => 'ncbi_sid=9FF573148F458940_CEC6SID; domain=.nih.gov; path=/; expires=Thu, 14 Mar 2019 22:36:03 GMT',
                                        'cache-control' => 'private',
                                        'access-control-allow-origin' => '*',
                                        'date' => 'Wed, 14 Mar 2018 22:36:03 GMT',
                                        'vary' => 'Accept-Encoding',
                                        'client-ssl-cert-issuer' => '/C=US/O=DigiCert Inc/OU=www.digicert.com/CN=DigiCert SHA2 High Assurance Server CA',
                                        'client-ssl-cipher' => 'ECDHE-RSA-AES256-GCM-SHA384',
                                        'client-peer' => '130.14.29.110:443',
                                        'strict-transport-security' => 'max-age=31536000; includeSubDomains; preload',
                                        'client-date' => 'Wed, 14 Mar 2018 22:36:03 GMT',
                                        'content-type' => 'text/xml; charset=UTF-8',
                                        'x-xss-protection' => '1; mode=block',
                                        'client-transfer-encoding' => [
                                                                        'chunked'
                                                                      ],
                                        'server' => 'Finatra',
                                        'client-ssl-socket-class' => 'IO::Socket::SSL',
                                        'client-response-num' => 1,
                                        'ncbi-phid' => 'DAF6EB0134DF5AEC039F11D248910BE7.1',
                                        'ncbi-sid' => '9FF573148F458940_CEC6SID',
                                        '::std_case' => {
                                                          'set-cookie' => 'Set-Cookie',
                                                          'access-control-allow-origin' => 'Access-Control-Allow-Origin',
                                                          'client-ssl-cert-issuer' => 'Client-SSL-Cert-Issuer',
                                                          'client-ssl-cipher' => 'Client-SSL-Cipher',
                                                          'client-peer' => 'Client-Peer',
                                                          'strict-transport-security' => 'Strict-Transport-Security',
                                                          'client-date' => 'Client-Date',
                                                          'x-xss-protection' => 'X-XSS-Protection',
                                                          'client-transfer-encoding' => 'Client-Transfer-Encoding',
                                                          'client-ssl-socket-class' => 'Client-SSL-Socket-Class',
                                                          'client-response-num' => 'Client-Response-Num',
                                                          'ncbi-phid' => 'NCBI-PHID',
                                                          'ncbi-sid' => 'NCBI-SID',
                                                          'x-ua-compatible' => 'X-UA-Compatible',
                                                          'client-ssl-cert-subject' => 'Client-SSL-Cert-Subject',
                                                          'content-security-policy' => 'Content-Security-Policy'
                                                        },
                                        'x-ua-compatible' => 'IE=Edge',
                                        'client-ssl-cert-subject' => '/C=US/ST=Maryland/L=Bethesda/O=National Library of Medicine/CN=*.ncbi.nlm.nih.gov',
                                        'content-security-policy' => 'upgrade-insecure-requests'
                                      }, 'HTTP::Headers' ),
                 '_msg' => 'Bad Request',
                 '_request' => bless( {
&rettype=full&retmode=xml&edirect_os=linux&edirect=7.90&tool=edirect&email=x@xxx.localdomain',
                                        '_uri' => bless( do{\(my $o = 'https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi')}, 'URI::https' ),
                                        '_headers' => bless( {
                                                               'user-agent' => 'libwww-perl/6.05',
                                                               'content-type' => 'application/x-www-form-urlencoded',
                                                               '::std_case' => {
                                                                                'if-ssl-cert-subject' => 'If-SSL-Cert-Subject'
                                                                               }
                                                             }, 'HTTP::Headers' ),
                                        '_method' => 'POST',
                                        '_uri_canonical' => $VAR1->{'_request'}{'_uri'}
                                      }, 'HTTP::Request' )
               }, 'HTTP::Response' );

What I've done:
Here is a sample of the data I used that worked correctly and yielded lineages when searching the taxIDs (taxIDs in 1st column):
549298,Contig1,Francisella halioticida,5.44e-126,gi|1213460577|gb|CP022132.1|,CP022132
240395,Contig2,Glandirana emeljanovi,4.74e-131,gi|1025812536|gb|KU641020.1|,KU641020
1069680,Contig3,Pneumocystis murina B123,2.34e-85,gi|1124120202|ref|XM_019757798.1|,XM_019757798
3871,Contig4,Lupinus angustifolius,1.89e-93,gi|1239725835|gb|CP023132.1|,CP023132

Here is a sample of the data that I am trying to get to work (when I search these manually against Taxonomy from their web page they seem to work and give lineages; it's just that I have just enough of them to make a non-programmatic solution really time consuming and tedious):
258234
1316279
1316315
710996
711000
225338
985895
710996
1382731
77133
344338

This is the business end of the script that I used that worked for the  working data:
#!/bin/bash
while IFS=',' read -r f1 f2 f3 f4 f5 f6
do
        #Firstly echo the Contig# (Query Seq. ID) followed by lineage information.
        echo -n "$f2," >> $2

        #Now get the lineage, pipe to xtract and parse XML for readability in output file
        efetch -db taxonomy \
        -id $f1 \
        -format xml \
        | xtract -pattern Taxon \
                -tab ',' -sep ',' \
                -def "NONE" \
                -element TaxId,ScientificName \
                -division LineageEx \
                -group Taxon \
                -if Rank -equals superkingdom \
                -or Rank -equals kingdom \
                ...
                -or Rank -equals species \
                -tab ',' -sep ',' \
                -def "NOT_AVAILABLE" \
                -element Rank,ScientificName \
                >> $2
done < "$1"

This is the slightly altered script that I am trying to use to make the not-working set of data yield results (and that produces the error):
#!/bin/bash
while IFS='' read -r line || [[ -n "$line" ]]; do
        efetch -db taxonomy \
        -id $line \
        -format xml \
        | xtract -pattern Taxon \
                -tab ',' -sep ',' \
                -def "NONE" \
                -element TaxId,ScientificName \
                -division LineageEx \
                -group Taxon \
                -if Rank -equals superkingdom \
                -or Rank -equals kingdom \
                ...
                -or Rank -equals species \
                -or Rank -equals subspecies \
                -tab ',' -sep ',' \
                -def "NOT_AVAILABLE" \
                -element Rank,ScientificName \
                >> $2
done < "$1"

As you can see, the only difference is that I am reading in entire lines for the TaxIDs list because it is simply a list with line breaks between the IDs as opposed to a true CSV with other values as in the example that works.  I've saved the taxIDs list as both a CSV and TXT and it doesn't make a difference.  Also: I can echo what $line is in the second script and it is exactly I'd expect--just the numbers of the taxID, so why would I be getting a 400 error for bad syntax on my query if it looks exactly like what I am sending in the first, working script?
Thanks!

Hello, nice question. Just, could you please tame ([edit]) a bit the error log? It's very hard to look look at.

Aside from your specific question slightly, you might be interested in looking at the NCBITaxa() tool in the ETE3 package. It basically does what your script does, but provides lots of other useful methods for querying the ncbi taxa database.

If I understand correctly, your question boils down to: parsing a bunch of numbers in a row from a .csv file gives a different result from parsing a bunch of numbers as rows of a text file.  I that right?  I suggest that you look at the source of your text file.  For example, was the text file created on a Windows host (running bash, or in the command terminal?)?
Try using this modified version of the unix cat command:
cat -et text_file.txt
The very end of each line will have a '\$' sign and that will show you if there are any extra non-visible characters in your rows.  For example, a DOS/Windows formatted file will usually have a ^M at the end of each row before the '$'.  If you don't strip off those line endings then they get added to your queries, even if you cannot see them.  You can strip off those linefeed (or are they carriage returns?) using a Perl one-liner like:
perl -pe 's/\r//' badfile.txt > goodfile.txt

Thanks, I suspect that between generating these files on my school's servers through bash to downloading them locally and manipulating them in Excel, that there must've been '$' characters like you said.  Since then, I've learned to do what I was doing in Excel on the command line and think I'll be able to avoid this altogether.

all_licenses: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Can I use DRIMseq to model interaction effect of variables?
I am trying to use DRIMseq for DTU with 2 treatments on two different strains of mice. In my previous analysis with genes only, I had found a strong interaction effect  of strain and treatments. Is it possible to model that with DRIMseq?
Also when I try to add strain term and use it as the batch effect as described in the tutorial, I get the following error:
! Using a subset of 0.1 genes to estimate common precision !
Error in optimHess(par = par, fn = dm_lik_regG, gr = dm_score_regG, x = x,  :
  non-finite value supplied by optim
In addition: There were 50 or more warnings (use warnings() to see the first 50)

Please see the sessionInfo() below.
R version 3.4.0 (2017-04-21)
Platform: x86_64-redhat-linux-gnu (64-bit)
Running under: CentOS Linux 7 (Core)

Matrix products: default
BLAS/LAPACK: /usr/lib64/R/lib/libRblas.so

locale:
 [1] LC_CTYPE=en_GB.UTF-8       LC_NUMERIC=C               LC_TIME=en_US.UTF-8        LC_COLLATE=en_GB.UTF-8     LC_MONETARY=en_US.UTF-8
 [6] LC_MESSAGES=en_GB.UTF-8    LC_PAPER=en_US.UTF-8       LC_NAME=C                  LC_ADDRESS=C               LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C

attached base packages:
[1] stats4    parallel  stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
 [1] DRIMSeq_1.6.0                            pheatmap_1.0.8                           ggplot2_2.2.1
 [4] edgeR_3.16.5                             limma_3.30.13                            TxDb.Mmusculus.UCSC.mm10.knownGene_3.4.0
 [7] GenomicFeatures_1.26.4                   AnnotationDbi_1.36.2                     Biobase_2.34.0
[10] GenomicRanges_1.26.4                     GenomeInfoDb_1.10.3                      IRanges_2.8.2
[13] S4Vectors_0.12.2                         BiocGenerics_0.20.0                      tximport_1.2.0

loaded via a namespace (and not attached):
 [1] Rcpp_0.12.12               RColorBrewer_1.1-2         plyr_1.8.4                 pillar_1.2.1               compiler_3.4.0
 [6] XVector_0.14.1             bitops_1.0-6               tools_3.4.0                zlibbioc_1.20.0            biomaRt_2.30.0
[11] digest_0.6.12              bit_1.1-12                 gtable_0.2.0               RSQLite_2.0                memoise_1.1.0
[16] tibble_1.4.2               lattice_0.20-35            pkgconfig_2.0.1            rlang_0.2.0                Matrix_1.2-10
[21] DBI_0.7                    stringr_1.3.0              rtracklayer_1.34.2         Biostrings_2.42.1          locfit_1.5-9.1
[26] bit64_0.9-7                grid_3.4.0                 XML_3.98-1.10              BiocParallel_1.8.2         magrittr_1.5
[31] reshape2_1.4.3             blob_1.1.0                 MASS_7.3-47                scales_0.5.0               Rsamtools_1.26.2
[36] GenomicAlignments_1.10.1   SummarizedExperiment_1.4.0 colorspace_1.3-2           labeling_0.3               stringi_1.1.7
[41] lazyeval_0.2.1             munsell_0.4.3              RCurl_1.95-4.10

One way to see the interaction effect would be including it in the design matrix as explained in the limma tutorial. Something on the line of:
design_full <- model.matrix(~ treatment*strain, data = samples(d))

Or you could build it yourself in a matrix like
 (Intercept) strainA_treatmentA strainA_treatmentB  strainB_treatmentA strainB_treatmentB
       1                      0                  1                  0                   0
       1                      1                  0                  0                   0

About the error, it seems like the data is not in the right format, you have either NA or Inf somewhere in the matrix. Check it and rethink the processing you made.


text:
How to verify if the stationary states of a biological model are correct?
I'm mathematically modeling yellow fever, and yet, after assembling the model, I need to find the stationary states of the model. I had to do it by hand because wolfram software could not get the necessary values.
My question is, how do I check if my stationary states are in line with my biological model?
I was told that there is a way to verify this by performing a substitution of the s values of the points within my own template code. But I have no idea how to do it.
The code used is taking forever to get back to me as this data behaves in my model

function ypsir = ypsirtry(t,y)
betah = .01;
muh = 0.01;
sigmah = 0.05;
alphah = 0.06;
deltah = 0.05;
mui = 0.07;
mum = 0.07;
betai = 0.03;
lambdai = 0.06;
lambdah = 0.5;
k = 500;

ypsir = zeros(5,1);
ypsir(1) = sigmah*y(3) - muh*y(1) - betah*y(1)*y(5) + lambdah*(y(2)+ y(3) +
y(1));
ypsir(2) = -deltah*y(2) - y(2)*(alphah + muh) + betah*y(5)*y(1);
ypsir(3) = deltah*y(2) - sigmah*y(3);
ypsir(4) = -mui*y(4) - mum*y(4) - betai*y(2)*y(4) + lambdai*(y(4) + y(5))*(1
- (y(4) + y(5))/k);
ypsir(5) = betai*y(2)*y(4) - mui*y(5) - mum*y(5);
end

Along with this calling program

clear;
to = 0;
tf =50;
yo = [99 1 0 1 3];
[t y] = ode45(@(t,y)ypsirtry(t,y),[to tf],yo);
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5))
title('First Model of Yellow Fever')
xlabel('time')
ylabel('susceptible, infected, recovered,Mosquito Non-carrier, Carrier
mosquito')

My stationary states are:

S = (1/betai)*[(muh+alphah)+(sigmah/deltah)*deltah]*(mui+mum)/betai*N*I);

I = sigmah*R/deltah;

R = ((muh/betah + lambdah)*(muh+alphah))/(sigmah+lamdah)-(muh+alphah)*
(deltah/sigmah*deltah)-(betah/muh - lambdah)*(betai*N*I)*
(delta/sigmah*deltah));

N = -k/lambdai((-lamdai+(lambdai/k))*P/-betai*I + lambdai - (2k/lambdai*P)-
mui-mum);

P = betai*N*I/(mui+mum);

At where
S = susceptible
I = infected
R = retweeted
N = non-carrier mosquito
P = carrier mosquito

I think your question is a good question, but we would need more details in order to answer them. How did you code your model? Which stationary states? (Why do you use the biomart tag?) Which data do you have to verify your models (if it is not purely theoretical)? What do you mean by "performing a schedule" ? You can [edit] it to include this or other information.

I'm already going to make the issues, sorry for the delay! I collared biomat because I thought it would be mathematical modeling in biology. I'm new to the site and I was a little lost on which site to put. I will post the code with the respective data. But when putting it to compile it is dying a lot of time. But I will add other information. Do you know how I can put on latex around here? to make it easier to see formulas in my edits. Sorry English, it is not my mother tongue how can you perceive rs ...

Wouldn't it be sufficient to determine this by running it until the system elicits stationary behavior and then comparing the values? If you did that for a variety of values for the constants (betah and so on) then that would also be a good test.

@wgrt See the [tour] page to learn about the site and the [help] page to know how it works. How much is a lot of time: 3 days? Do you know if it can be optimized? (I don't know enough matlab to propose it myself). Also you left some questions from my previous comment unanswered

@Llopis You're not going to optimize an RK45 method that's implemented in C much. If this really is slow to compute then presumably there are some very drastic changes in values (RK45 shrinks its steps in such cases).

@DevonRyan sometimes the methods are optimized but one doesn't use them efficiently. Anyway, wgrt: any update?

After much analysis of the model, I came to the conclusion that I will not be able to verify how the equilibrium points behave graphically. Since you will have to create a phase diagram in R ^ 5, there are five equations, and I can not do two-to-two comparison in these cases, because I will always have a function that depends on a variable that will not be placed, but will be in the equation plotted anyway.
But in any case, if it were in three dimensions, for the purpose of information for those who are interested, the verification of the points would be to place the points (which in calculating have both numbers and not functions in function of other functions), within the model created, hence plot the phase diagram and check whether the arrows converge or not to the point. If they are converging, they are correct.
Regarding the model is taking forever, I took the logistics part and ran well!


text:
Digisome generation
I scoured and searched the web for some implementation of digisome creation, but I had a hard time finding any good source code.
By digisomes, what I mean is images with chromosomal masks/borders that have been simulated. I am seeking this because I am trying to do a pixel level prediction on chromosomes using a deep network. And, I need ground truths with the chromosomes being illuminated. So ideally, I would like a binary image with just the chromosomes being illuminated.
How exactly are digisomes generated? Is a single actual chromosome extracted from some metaphase image, then morphed in interesting ways, and finally, all the different morphs are stitched together in one big collage format so that it there are numerous chromosomes all seemingly different?
Also, if someone could point me to some good source code, that'd be helpful.

Please [edit] your question and explain what digisomes are. I found one reference to them from the proceedings of a conference from 1990 and they seem to be referring to some form of in silico simulation of chromosomes. Is that what you mean? If the term is in use, it is not very common outside whatever specific subfield the authors may have been working in. Also, please explain what sort f analysis you would want to use these for since the kind of simulation/model you will need depends entirely on the questions you want to answer.

I've edited my question. Please let me know if that's not sufficient.

That's much better, thanks. I still don't get what you need, but now I'm sure that's just because of my own ignorance of the field. I also edited so that the extra information you added is integrated into the question. Since there's an edit history here, and the questions can (we hope) stick around for years, it's better to edit so that the question makes sense to someone reading it for the first time.

I found a very nice github repo that simulates chromosomes and their labels for overlapping chromosomes - and with a little bit of tweaking, it can be amended for a variety of different simulations.
This was also used in a Kaggle post; so right up the alley for deep learning.


text:
How are positions visually indicated in the UCSC Genome Browser?
I'm a little confused about how the position indices work in the browser. Based on the picture below, is base number 755440 an A or a G?
The general question is how do the indexing guides work? Is the base right under the | the one corresponding to the index or is it the one just before the |?

It's an A. The number refers to the position just before the |. This is clearer on different zoom levels, but I double checked by going to the very beginning of a chromosome:

As you can see in the image above, the | is always after the base it refers to. I also confirmed by checking a few random positions using another genome browser which places the positions unambiguously above each base.
For more details, see this UCSC blog post.

Yep and the reason for this is because we store everything internally as zero based coordinates but then display everything in a 1-based format. For more information see the following blog post: http://genome.ucsc.edu/blog/the-ucsc-genome-browser-coordinate-counting-systems/

@ChristopherLee thank you! I had actually written out an answer stating that it must be the position after since bed files start at 0 and then thought to check, just to make sure, and was very much surprised that it was actually the position before!


text:
Pepsin digest (cleavage) does not work using RE?
Aim is to code the theoretical peptic cleavage of protein sequences in Python. The cleavage (cutting) rule for pepsin is: 1234^567. This is a position, ^ stands for the cleavage point.
There are rules that it mustn't have at all and rules where it must have at least one. The rule for mustn't is in position 2, it mustn't be R, K or H. In position 3 and 7, there can't be a P and in position 4, there mustn't be an R.
If it follows "the must not" rules, then it will check "the must" rules so that if any of one them is met, it will cleave. So if there is a P in position 1 or 2, or if there's an L, W, Y or F in position 4 or 5 then, it will be cleaved (use newline).
This is the code I have so far, but it doesn't work yet:
import re
print('INSERT YOUR PEPTIDE')
X= input();
r=re.findall('(?<=[P])(?<=[P,^R,^K,^H])(?<=([^P])(?<=[L,W,F,Y,^R])(?<=[L,W,Y,F])/D(?<=[^P])',X)
print(r)

My input is:
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

How can I make this work?

I tried to clarify your question, please check and make sure I didn't change anything (it was very hard to understand). It would be much easier for us to help you if you could give some examples of input matching different rules and the output you expect from each of them.

I am really not sure if I've understood your rules correctly, but I think you want this:

If all of these are true:

position 2: not R, K or H
position 3: not P
position 4: not R
position 7: not P

And, if at least one of these is true:

position 1: P
position 1: P
position 4: L or W or Y or F
position 5: L or W or Y or F

Then cleave.
If so, then you probably don't want to do it using regular expressions. First because it will be hard, and second because it will be absolutely awful to maintain if you need to come back to it. Instead, try something like this:
#!/usr/bin/env python3
# X = input();
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
seq = list(X)
for i, res in enumerate(seq):
    if i > len(seq) - 7:
        break
    if (
            seq[i+1] in ('R','K','H') or
            seq[i+2] == 'P' or
            seq[i+3] == 'R' or
            seq[i+6] == 'P'
    ):
        continue
    if res == 'P':
        start = i
    elif seq[i+1] == 'P':
        start = i+1
    else:
        continue
    if (
            seq[i+4] in ('L', 'W', 'Y', 'F') or
            seq[i+5] in ('L', 'W', 'Y', 'F')
    ):
        print("%s\n%s\n%s" %
            ("".join(seq[0:start]),
             "".join(seq[start:start+7]),
             "".join(seq[start+7:])
            )
        )

This will produce the following output when run:
$ foo.py
TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMV
PKCFNEL
LMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI

If you insist on using regular expressions, you could separate it into two expressions and use something like this:
import re
X = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
re1='P[^RKH][^P][^R][LWYF].[^P]'
re2='P[^P][^R][LWYF].[^P]'
r=re.findall('(?:P[^RKH][^P][^R][LWYF].[^P])|(?:P[^P][^R][LWYF].[^P])',X)
print(r)

Before you do, however, please consider the wise words of XKCD 208:

thank for your reply but Okay I confuse with the rule so it should be like in this follow 2 pattern
pattern1
position 2: not R, K or H
position 3: not P
position 4: not R
position 5: F or L or W or Y
position 6: not P
pattern2
position 2: not R, K or H
position 3: not P
position 4: F or L or W or Y
position 6: not P

I just made the pattern1 but still no idea how to combine two of these pattern
import re
myprotein = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
my_protein = []
for protein in myprotein:
myprotein = protein.rstrip('\n')
my_protein.append(myprotein)
my_pro = (''.join(my_protein))
#cleaves sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])','\n', my_pro)
print (peptides)
It's must be great if you could help.

Ok I think I almost solve this problem so I still using RE
import re
myprotein = "TTGWMFKGLYQLKSDEARMHNLQRWYVFKVIEFWGINRKDQHDSRYWFSHFCWKLEHGMVPKCFNELLMPPMSTHRSMFNAYSMSFCVQVKEWRGSISRMCDHHQQQCQYWRLMAGGTLI"
my_protein = []
for protein in myprotein:
    myprotein = protein.rstrip('\n')
    my_protein.append(myprotein)
    my_pro = (''.join(my_protein))

   #cleaves sequence
peptides = re.sub(r'(?<=[^HRK][^P][^R])(?=[FLWY][^P])|(?<=[^HRK][^P][FLWY])(?=.[^P])','\n', my_pro)

print (peptides)

and it produce to
TTG
W
M
F
KG
LY
Q
L
KSDEARMHN
LQRW
YV
F
KVIE
F
W
GINRKDQHDSRY
WF
SH
F
CW
K
L
EHGMVPKC
FNE
L
LMPPMSTHRSMF
NA
Y
SMS
F
CVQVKE
WRGSISRMCDHHQQQCQ
Y
W
RL
MAGGT
LI

but the true answer should be like (checking from https://web.expasy.org/peptide_cutter/ )
TTG
W
M
F
KG
LY
Q
L
KSDEARMHN
LQRW
YV
F
KVIE
F
W
GINRKDQHDSRY
WF
SH
F
CW
K
L
EHGMVPKC
FNE
L
LMPPMSTHRSMF
NA
Y
SMS
F
CVQVKE
WRGSISRMCDHHQQQCQ
Y
W
RL
MAGGT
L
I

so how can I make this work?

In addition if I manually add some space after the sequence  it will get the right result.

You don't need the space, just add a * to the last [^P]: [^P]* so it matches "0 or 1 not P". But don't you need a P at position 1 or 2?

Oh thanks! I found that it is not necessary for P at position 1 or 2 so it's okay for no  P at position 1 or 2 :)

Ah, sorry, make that a ? for "0 or 1", instead of a * for "0 or more"

thank you kindly oh I have a few more question. so if I want to make it to loop that cleavage more than one sequence till the last sequence and every sequence start with ">" above (fasta format) and I want result that separate each sequence how can I make it work.

fasta format example:

sequence1
ATGCGGAGCTTAGATTCTCGAGATCTCGATATCGCGCTTATAAAAGGCCCGGATTAGGGC
TAGCTAGATATCGCGATAGCTAGGGATATCGAGATGCGATACG
sequence2
GTACTCGATACGCTACGCGATATTGCGCGATACGCATAGCTAACGATCGACTAGTGATGC
ATAGAGCTAGATCAGCTACGATAGCATCGATCGACTACGATCAGCATCAC


text:
Hierarchial PCA Clustering with duplicated row names
I have a matrix object in R of 1500 rows and 20 columns. Based on my own algorithm of analysis, I know there are roughly 7 clusters. In my matrix object, each row is one data point, and is labeled to which cluster it belongs by its row name. Thus, there are duplicated row names. I wanted to see if PCA would give me the approximate same clustering of my data.
Segment1 Segment2 Segment3 Segment4 Segment5 Segment6 Segment7 Segment8 Segment9 Segment10 Segment11 Segment12 Segment13 Segment14 Segment15 Segment16 Segment17 Segment18 Segment19 Segment20
B  38.40000 41.75200 44.38400 44.18400 45.37600 37.49600 41.36800 33.93600 38.00800  42.51200 46.49600  40.48000  45.40800  46.32800  43.78400  39.88800  38.84000  40.56800  42.03200 38.89185
C  45.53846 50.08462 39.91538 36.95385 34.96154 39.74615 38.01538 35.75385 35.54615  36.69231 35.20769  38.05385  39.29231  37.96923  37.30000  36.86923  39.19231  38.81538  43.69231 38.06400
C  46.05176 41.69412 38.80000 37.75529 39.67529 39.07765 39.17647 38.24941 39.58588  38.63529 38.30588  41.87765  38.97412  40.13647  42.27294  38.24471  35.41647  40.80000  38.07059 42.11294
B  44.20000 43.42857 44.80000 35.20000 35.91429 37.82857 51.45714 44.68571 46.68571  48.74286 41.25091  39.45455  38.17091  40.70182  40.39273  41.28727  40.63636  41.50909  41.68364 41.29455
E1 45.06909 41.09818 40.02909 42.50182 42.34909 39.84727 41.42909 40.47273 40.28000  40.51636 41.25091  39.45455  38.17091  40.70182  40.39273  41.28727  40.63636  41.50909  41.68364 41.29455
E3 40.87407 39.27704 44.13630 43.25037 35.86667 37.30667 38.76148 40.74667 38.93333  43.16148 37.47259  37.73630  38.34370  39.00148  36.96889  37.76593  39.14667  37.92593  37.62963 38.89185

(The first column are the row names, either being A,B,C,D,E1,E2 or E3 according to my clustering)
I used various packages and functions like pcomp(), princomp() from the stats package, and PCA() from the FactoMineR package. If I wanted to do hierarchical clustering with the HCPC() function from the same package, it was not possible because of the duplicated row names.
Using:
fviz_pca_ind(pca1, geom.ind = "point", axes = c(1,3) , col.ind = new$Class, palette = c("#66C2A5" ,"#FC8D62" ,"#8DA0CB", "#E78AC3" ,"#A6D854", "#FFD92F", "#E5C494"), addEllipses = T, legend.title="Classes")

I got the individual datapoints on my two PC axes with ellipses around them.
But for clustering:
PCA_FR <- PCA(mat_data, ncp = 4, graph = FALSE)
Hier_PCA<- HCPC(PCA_FR, graph = F)

I couldn't compute Hier_PCA because of duplicated row names.
> Hier_PCA<- HCPC(PCA_FR, graph = F)
Error in `contrasts<-`(`*tmp*`, value = contr.funs[1 + isOF[nn]]) :
contrasts can be applied only to factors with 2 or more levels
In addition: Warning messages:
1: In data.row.names(row.names, rowsi, i) :
some row.names duplicated:  3,4,7,8,10,11,12,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,

How can I best cluster my data using the first three PCs, and visualize the 1500 datapoints, but labeled only in 7 ways?

Did you tried to use unique rownames and store the name of the cluster each data belongs in another vector? What was then the problem?

I didn't try that because I want to see my clusters according to the labels A,B,C,D,E1,E2,E3. Doing that, I thought I would still see just numbers, and not my labels.

Yes, but once the position are calculated you can use the base function text to plot the labels at that positions as previously commented. (Please do not remove the comments so quickly moderator :)

Aah okay. Can you show me the code on how to use the text function and separate my rows according to the rownames, in the answer?

I could in a while, but search in StackOverflow, I'm sure this has been answered somewhere (I  learned from there)

Try to make your rownames unique first.
So if you have this matrix with duplicate rownames:
m <- replicate(20, rnorm(6))
rownames(m) <- c("B", "C", "C", "B", "E1", "E3")
m
         [,1]        [,2]      [,3]      [,4]       [,5]        [,6]       [,7]
B   2.4071997  1.30040528 1.0765500 0.2005440  1.1732384  2.34462132  0.3702172
C  -0.2683534 -0.74607893 0.1611073 0.3145887  0.8208741  1.26575765  0.4205310
C   0.1997264  0.03318396 0.1967151 0.7785350  0.9382146 -0.32475195 -0.9832998
B  -1.2726538  0.22519992 1.1556121 0.7161238  0.5601368  1.02900724 -0.2726709
E1  0.8062722 -1.80899128 0.3936155 0.5724118 -0.8163454  0.04306317 -0.4074828
E3  0.3320706  0.21218846 0.9210126 0.9151837  0.4588621 -1.27784583  1.3292399
         [,8]       [,9]       [,10]      [,11]       [,12]       [,13]
B   0.2305559  0.2964808 -0.01035426 -2.6062776 -0.03143381 -1.48828554
C   0.1622548 -0.6244459 -0.66571655 -0.3529994  1.07540375 -1.58949697
C   0.5984752 -0.6423845 -1.00001700  0.2212160  0.05499773  1.36182121
B  -0.5594747 -0.3570748  1.44098577 -0.4939490 -0.93386684  0.54919399
E1  1.6378047  0.8702920  3.06641490  0.1948452  0.36339348 -1.04694826
E3 -2.2422997  1.0682864  0.14326146 -0.4513848 -0.25819231  0.01364221
        [,14]      [,15]      [,16]      [,17]       [,18]      [,19]
B  -2.0410367  2.4858879  1.0564085 -1.7353230 -0.63024160 -1.2642884
C  -1.2216387  0.5533529 -0.2757584 -0.2119105 -0.68899706  1.0587382
C   0.6401224 -1.5998479  1.4072007 -1.5062900 -0.60161624  0.9781313
B   0.4309208  1.4310142 -0.1702066 -0.7195977  0.37466185  0.2605195
E1  1.2538412  0.3413762 -1.9346107  0.7314531 -0.03105008  0.7341104
E3 -0.3146694 -1.3028086  0.5043235 -0.2568268  0.08336845  1.3097557
        [,20]
B  -1.1245496
C  -0.4552617
C  -1.5498408
B   0.1261818
E1 -0.1180164
E3  2.1035368

# Make them unique with make.names function

rownames(m) <- make.names(rownames(m), unique=T)

m
          [,1]        [,2]      [,3]      [,4]       [,5]        [,6]
B    2.4071997  1.30040528 1.0765500 0.2005440  1.1732384  2.34462132
C   -0.2683534 -0.74607893 0.1611073 0.3145887  0.8208741  1.26575765
C.1  0.1997264  0.03318396 0.1967151 0.7785350  0.9382146 -0.32475195
B.1 -1.2726538  0.22519992 1.1556121 0.7161238  0.5601368  1.02900724
E1   0.8062722 -1.80899128 0.3936155 0.5724118 -0.8163454  0.04306317
E3   0.3320706  0.21218846 0.9210126 0.9151837  0.4588621 -1.27784583
          [,7]       [,8]       [,9]       [,10]      [,11]       [,12]
B    0.3702172  0.2305559  0.2964808 -0.01035426 -2.6062776 -0.03143381
C    0.4205310  0.1622548 -0.6244459 -0.66571655 -0.3529994  1.07540375
C.1 -0.9832998  0.5984752 -0.6423845 -1.00001700  0.2212160  0.05499773
B.1 -0.2726709 -0.5594747 -0.3570748  1.44098577 -0.4939490 -0.93386684
E1  -0.4074828  1.6378047  0.8702920  3.06641490  0.1948452  0.36339348
E3   1.3292399 -2.2422997  1.0682864  0.14326146 -0.4513848 -0.25819231
          [,13]      [,14]      [,15]      [,16]      [,17]       [,18]
B   -1.48828554 -2.0410367  2.4858879  1.0564085 -1.7353230 -0.63024160
C   -1.58949697 -1.2216387  0.5533529 -0.2757584 -0.2119105 -0.68899706
C.1  1.36182121  0.6401224 -1.5998479  1.4072007 -1.5062900 -0.60161624
B.1  0.54919399  0.4309208  1.4310142 -0.1702066 -0.7195977  0.37466185
E1  -1.04694826  1.2538412  0.3413762 -1.9346107  0.7314531 -0.03105008
E3   0.01364221 -0.3146694 -1.3028086  0.5043235 -0.2568268  0.08336845
         [,19]      [,20]
B   -1.2642884 -1.1245496
C    1.0587382 -0.4552617
C.1  0.9781313 -1.5498408
B.1  0.2605195  0.1261818
E1   0.7341104 -0.1180164
E3   1.3097557  2.1035368


text:
Extracting strand specific reads from MinION cDNA-PCR protocol
I recently performed my first MinION run, and now I'm trying to analyze the data. Being pretty new to the bioinformatics field, I was hoping some of you could help me out.
As a bit of background, I'm trying to define transcript isoforms, and down the line, I'm hoping to quantify isoform abundances. As my first attempt, I used a 1D kit that first performs RT starting at the polyA tail and then template switching to enrich for full length transcripts. To initiate RT, a poly-dT primer is added and after template switching a second primer is added. Using those added primer sequences, I then performed PCR to amplify the full-length cDNAs. This seemed to work well, so I went ahead and sequenced the library.
Now, on to my problem. I have aligned my fastq file to my genome with both bwa-mem and minimap2. In both cases, I get alignment of reads to the genome where I expect transcription to be occurring, but half of the reads are in the sense direction and half are in the antisense direction. This is presumably due to the fact that the cDNA is double stranded and either strand can be recruited to pore to be sequenced. However, because I'm interested in mapping transcript isoforms which could include anti-sense transcripts, this isn't exactly what I hoped for.
I'd like to figure out a way to extract only the reads that correspond to the actual mRNA transcript's sequence and not the reverse complement. Does anyone have ideas as far as programs that would help with this?
I'm thinking that I should be able to extract out only the molecules that have the polyA or have the strand switching primer sequence or something like that, but I'm not sure how to properly implement this.

This has been set up as a bioinformatics protocol on protocols.io.
I've been using LAST to identify the adapter orientation relative to the genome, and then using that information to split BAM files up and recombine them to create two strand-specific files that are displayable in a genome browser.
I have written my own script to process LAST results into a CSV format, which makes it easier to do line-by-line data filtering.
Note: these scripts have been slightly modified from scripts that I have run. Consider them demonstrative: pay attention to the words rather than the script.
Read correction
I prefer starting off my data analysis with a read correction with Canu (ideally v1.7 now that it's out, because that attempts correction of all reads, but here I use canu v1.6). I use minimap as the mapper to speed this up. The genomeSize parameter should be approximately a tenth to a fortieth of the number of bases in your dataset to make sure that no sequences are excluded (bigger is better, as long as Canu doesn't freak out about memory consumption):
~/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

This creates a file 4T1_BC06/4T1_BC06.correctedReads.fasta.gz.
Chimeric read filtering
The next step I carry out is a basic read-level QC to exclude chimeric reads. Porechop can be used for this, although that removes adapters by default, which is not particularly useful in this case.
I use LAST to search for adapter sequences within the corrected reads, pass it through my conversion script, and extract out duplicated mappings (i.e. where the same read/adapter pair appears more than once in the mapping results):
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06/4T1_BC06.correctedReads.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | awk -F',' '{print $1,$2}' | sort | \
  uniq -d | awk '{print $1}' | uniq > reads_with_duplicated_adapters.txt

I use this file to filter out chimeric reads from the corrected dataset using another fastq filtering script I've created:
~/scripts/fastx-fetch.pl -v -i reads_with_duplicated_adapters.txt 4T1_BC06/4T1_BC06.correctedReads.fasta.gz | \
  gzip > 4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Adapter filtering
Along roughly similar lines to the chimeric read filtering, I then look for the strand switch and VNT adapters in the sequences. For a forward-orientated query, I expect the strand switch primer to be in the forward direction, and the VNP primer to be in the reverse direction:
lastal -P 10 ONT_barcodes_adapters.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,+' -e 'ONT_VNP,-' | \
  awk -F',' '{print $1}' | sort | uniq > fwdQry_seqs_BC06.txt
lastal -P 10 barcodes_primerSeqs.fa <(zcat 4T1_BC06.correctedReads.uniqueOnly.fasta.gz) | \
  ~/scripts/maf_bcsplit.pl | grep -e 'ONT_SSP,-' -e 'ONT_VNP,+' | \
  awk -F',' '{print $1}' | sort | uniq > revQry_seqs_BC06.txt

And then more filtering to split the reads up into forward and reverse [genome-relative] subsets:
~/scripts/fastx-fetch.pl -i fwdQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz
~/scripts/fastx-fetch.pl -i revQry_seqs_BC06.txt 4T1_BC06.correctedReads.uniqueOnly.fasta.gz | \
  gzip > rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz

Mapping
Now that the reads have been oriented, the mapping can be done:
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  fwd_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > fwd_4T1_BC06.CU_vs_mmus.bam
~/install/minimap2/minimap2 -t 10 -a mmus_ucsc_all_cdna.idx -x splice \
  rev_4T1_BC06.correctedReads.uniqueOnly.fasta.gz \
  > rev_4T1_BC06.CU_vs_mmus.bam

Splitting and recombining
The mapped files are then split based on their mapping direction using samtools view flag filtering to exclude or include reverse-mapped reads:
samtools view -F 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > fwd_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b fwd_4T1_BC06.CU_vs_mmus.bam > rev_fwd_4T1_BC06.CU_vs_mmus.bam
samtools view -F 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools view -f 0x10 -b rev_4T1_BC06.CU_vs_mmus.bam > rev_rev_4T1_BC06.CU_vs_mmus.bam

And finally the files are recombined to identify forward and reverse-encoded transcripts via samtools merge. I'm calling the "pos" strand the one that is encoded in the same direction as the genome, and the "neg" strand the reverse-complement direction:
samtools merge pos_4T1_BC06.CU_vs_mmus.bam \
  rev_fwd_4T1_BC06.CU_vs_mmus.bam fwd_rev_4T1_BC06.CU_vs_mmus.bam
samtools merge neg_4T1_BC06.CU_vs_mmus.bam \
  fwd_fwd_4T1_BC06.CU_vs_mmus.bam rev_rev_4T1_BC06.CU_vs_mmus.bam

Alternatively, it's possible to set the BAM flags for read 1 and read 2 for the reverse and forward-encoded adapters respectively, which allows everything to then be combined into a single BAM file (and treated in the same way as a strand-specific Illumina BAM file).
If this has worked properly, then mapping human or mouse to the mitochondrial genome should show most expression appearing on the positive strand, with a small scattering of negative-strand expression, a bit like this:

For a bit of a background on how/why this works, see this answer.


text:
Merging regions according to their identifier
I have a huge file (20 GB) which has a range of genomic locations, and for each location there is an identifier(4th column), which is sometimes the same.
file1.txt

chr1  10  20  ABC
chr1  13  20  ABC
chr1  14  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  21  JJJ
chr2  21  31  JJJ
chr2  23  27  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  101 MMM
chr3  101 105 MMM

I am trying to merge the ranges of column 2 and 3, if they are in the same chromosome and having the same identifier ( having the 1st and 4th columns the same).
I have tried to use bedtools, as shown in this post, which looks very similar but not the same with mine, as they are looking for different outputs ( range vs merge )
So the answers in that post was with using groupby in bedtools2,source code from github.
I have applied the same as the following:
./groupBy -i ../file1.txt -g 1,4 -c 2,3 -o min,max | awk -v OFS='\t' '{print $1, $3, $4, $2}'

chr1    10      37      ABC
chr2    15      56      JJJ
chr2    25      26      MMM
chr3    32      105     MMM

But it only connects the ranges of identifiers.. It groups by, and does not merge.
And, if I use bedtools v2.26 merge function :
>sort -k1,1 -k2,2n file1.txt > file2.txt
>cat file2.txt
chr1    10      20      ABC
chr1    13      20      ABC
chr1    14      21      ABC
chr1    22      27      ABC
chr1    29      37      ABC
chr2    15      21      JJJ
chr2    21      31      JJJ
chr2    23      27      JJJ
chr2    25      26      MMM
chr2    35      56      JJJ
chr3    32      42      MMM
chr3    45      76      MMM
chr3    88      101     MMM
chr3    101     105     MMM

bedtools merge -i file2.txt
chr1    10      21
chr1    22      27
chr1    29      37
chr2    15      31
chr2    35      56
chr3    32      42
chr3    45      76
chr3    88      105

Which, merges without seeing the identifiers as it also gives me this output.
What I want is to merge the ranges if they overlap, if having the same 1st and 4th column as:
expected_output.txt

chr1  10  21  ABC
chr1  22  27  ABC
chr1  29  37  ABC
chr2  15  31  JJJ
chr2  35  56  JJJ
chr2  25  26  MMM
chr3  32  42  MMM
chr3  45  76  MMM
chr3  88  105 MMM

I only have access to 15GB of RAM, so I can't store the entire file in memory.

Possible duplicate of Merging bed records based on name

This post is already linked in my entry, and they are looking for different outputs,so not a duplicate

The following python script will do what you want and should be relatively memory efficient. It processes a single chromosome at a time, so either sort the input or ensure that at least entries in a group are next to each other.
#!/usr/bin/env python
import sys

# Process a sorted list of tuples (e.g., [(0, 10), (10, 20), ...])
# Overlapping items will be merged and a list of the merged regions is
# returned (e.g., [(0, 20), ...])
def merge(values):
    start, end = values[0]
    o = []
    for v in values[1:]:
        if v[0] <= end:
            if v[1] > end:
                end = v[1]
        else:
            o.append((start, end))
            start, end = v
    o.append((start, end))
    return o

# Process a single chromosome, one group at a time
def processChrom(d, chrom):
    for k, v in d.items():
        v.sort()
        tuples = merge(v)
        for t in tuples:
            print("{}\t{}\t{}\t{}".format(chrom, t[0], t[1], k))

def main():
    currentChrom = None
    d = dict()
    for line in sys.stdin:
        cols = line.split()
        if cols[0] == currentChrom:
            if cols[3] not in d:
                d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
        else:
            if currentChrom:
                processChrom(d, currentChrom)
            d = dict()
            d[cols[3]] = []
            d[cols[3]].append((int(cols[1]), int(cols[2])))
            currentChrom = cols[0]

    if currentChrom:
        processChrom(d, currentChrom)

if __name__ == "__main__":
    main()

Usage is python foo.py < foo.txt and the output is:
chr1    10  21  ABC
chr1    22  27  ABC
chr1    29  37  ABC
chr2    25  26  MMM
chr2    15  31  JJJ
chr2    35  56  JJJ
chr3    32  42  MMM
chr3    45  76  MMM
chr3    88  105 MMM

You'll notice that this isn't necessarily sorted, so you may want to pass the output through bedtools sort or sort.

Thank you for your answer. Do you mean that the input file should be sorted by 1st column or 4th?

You might as well sort it by the first column. You could alternatively rewrite the program to handle files sorted by the fourth column, but the current script isn't written to handle that.

Sure. The code is clean and working well. I am glad that it is adaptable

If I understand correctly, you could create dummy chromosomes made by merging chromosome and identifier, merge with bedtools, split back chromosomes and identifiers. E.g.
awk -v OFS='\t' '{print $1"_"$4, $2, $3}' file1.txt \
| mergeBed \
| sed 's/_/\t/' \
| awk -v OFS='\t' '{print $1, $3, $4, $2}'

(Assuming "_" is not part of chromosome or identifier name, adjust accordingly if so).

Thank you for the answer, it works well. However, would bedtools complain if the $4th column has special characters within?

@bapors why would the 4th column have special characters?

@bapors if you mean that bedtools would complain if the modified 1st field would have special characters, then please edit your question to give a representative example of your input data. That said I just tried using the characters &^%*()@#^ and bedtools had no problem.

No, my input does not have special characters, my question was solely theoretical. Perfect!

Just to add this can be done with bedtools by processing each identifier separately which may or may not be practical given your large file.
for ID in $(perl -ane 'print $F[3] . "\n";' file1.txt | sort | uniq);do
    grep $ID file1.txt | \
    perl -ane 'chomp;print join("\t",@F[0..$#F]) . "\n";' | \
    sort -V | \
    bedtools merge -i stdin -c 4 -o distinct
done | sort -V

Here is a fairly simple (and hopefully readable) native Python solution. It assumes that the bed file is sorted prior to parsing:
from itertools import groupby

class BedRecord(object):
    def __init__(self, chr, start, end, feature):
        self.chr = chr
        self.start = start
        self.end = end
        self.feature = feature

def parse_records(in_file):
    with open(in_file) as f:
        for line in f:
            yield BedRecord(*line.strip().split())

def merge_groups(group):
    output = [next(group)]
    for record in group:
        if record.start <= output[-1].end:
            output[-1].end = max(output[-1].end, record.end)
        else:
            output.append(record)
    return output

def main(in_file):
    for _, group in groupby(parse_records(in_file), lambda x: (x.chr, x.feature)):
        for r in merge_groups(group):
            print('\t'.join([r.chr, r.start, r.end, r.feature]))

if __name__ == '__main__':
    main('test.bed')

Thank you for your answer. It assumes that the input will be sorted by which column?

@bapors effectively by chromosome, then feature then start value, you could easily add code to sort based on these criteria if this assumption is violated


text:
How to read gene regulatory network edge list files?
The following is an excerpt from an edge list file from the Gene Regulatory Database, YeastNet v3
YML100W YMR261C 5.73968189451551
YDR074W YML100W 5.73143623470086
YDR074W YMR261C 5.72385204041689
YML106W YMR271C 5.64698433340441
YGR240C YMR205C 5.64063955646679
YDL131W YDL182W 5.63907349642605
YNL037C YOR136W 5.62782704009976
YHR208W YJR148W 5.61712512462848
YAL062W YOR375C 5.59956647235823
YNL104C YOR108W 5.59046546303091
YGR132C YGR231C 5.58996626839351
YBR140C YOL081W 5.58738025271142
YDR453C YML028W 5.58281825879471
YNL096C YOR096W 5.58060076040965
YGR180C YJL026W 5.57785446902287
YGR056W YLR357W 5.57582611528462
YBR009C YNL030W 5.57261292043205
YMR169C YMR170C 5.55866364259678
YAR071W YBR093C 5.55668852247789
YER081W YIL074C 5.55252000887001

What do these rows mean?

It also depends which file you have downloaded. Apparently the evidence can be based on various inputs such as co-citation in PubMed, or co-occurrence of domains... The details are listed on their webpage.

You can read about it on the website, the last column is the Log likelihood scores (LLS) between the genes in the first two columns. An LLS of 0 means not better than random link.
On the website they show papers with fuller details about the statistics.


text:
Remove/delete sequences by ID from multifasta
I have a fasta file like this:
>Id1
ATCCTT
>Id2
ATTTTCCC
>Id3
TTTCCCCAAAA
>Id4
CCCTTTAAA

I want to delete sequences that have the following IDs.
Id2
Id3

The IDs are in a .txt file, and the text file will be used to match and delete those sequences.
My output should be a fasta file like this
>Id1
ATCCTT
>Id4
CCCTTTAAA

But I want it with awk and/or sed and/or bash (no python or perl).
How can I do it?

why the language restrictions?

Exactly: use the right tool for the job. Unless you can make some very strong assumptions about your input files, use a proper FASTA parser, not some hacked-together cruft.

Use whatever tool you want to use if it gets your work done before a deadline. A lot of parsers are slow as mud as compared with Unix tools, but I admit they take a lot of drudgery out of understanding how things work.

Suppose you keep sequence names in ids.txt and sequences in seq.fa:
awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa

Just today I wrote a script to do exactly this using Biopython. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script filter_fasta_by_list_of_headers.py is :
#!/usr/bin/env python3

from Bio import SeqIO
import sys

ffile = SeqIO.parse(sys.argv[1], "fasta")
header_set = set(line.strip() for line in open(sys.argv[2]))

for seq_record in ffile:
    try:
        header_set.remove(seq_record.name)
    except KeyError:
        print(seq_record.format("fasta"))
        continue
if len(header_set) != 0:
    print(len(header_set),'of the headers from list were not identified in the input fasta file.', file=sys.stderr)

and usage :
filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter > filtered.fasta

P.S. it's quite easy to turn over the script to extract the sequences from the list (just the print line would have to move after line header_set.remove(seq_record.name)
I made a modification to use in and benchmarked it on normal insect genome assembly - 1Gbp, 300,000 scaffolds and I gave it lists of 5,000 and 100,000 scaffolds to filter. Regardless of the size of the list or algorithm used it took 29-30s. In other words, if LBYL is faster than EAFP, it's on a scale that totally does not matter.

@KamilSJaron Nice!

Ah +1 for showing my intuitions to be wrong, one other thought, it should be faster to write all the records at once rather than one by one, is something like this faster? https://pastebin.com/n7RUgW5x Sorry if it seems like micro-optimization, I'm just curious

If you want to learn how to do things with command-line tools, you can linearize the FASTA with awk, pipe to grep to filter for items of interest named in patterns.txt, then pipe to tr to delinearize:
$ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" > out.fa

This will work on multiline FASTA.

Excellent solution, except it will keep the sequence specified in patterns.txt, not remove them

Good point. I must have missed that part of the question. One can use grep -v to remove elements that match specified patterns.

Using the FastaToTbl and TblToFasta scripts I have posted before, you can do:
FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa

Presume grep is OK if sed, awk are. The -A(n) and -B(n) flags give (n) lines post match. Presuming all your fasta sequence will be one line is a bit dangerous, but it works for your example. First get those IDs to be removed (in rmid.txt), then inverse grep against the initial fasta.
grep -A1 -f rmid.txt fasta.fa > rmfile.fa
grep -v -f rmfile.fa fasta.fa

The real answer is to use a script that defines a delimiter other than newline "\n", then parse for IDs, so better to use one of the languages you don't want to use...

This assumes only one line of sequence per ID which is not likely true for fasta files (but is often true of fastq).

The problem with greping Id1 is substrings, this will also capture Id10 etc. you could add --word-regexp flag, but this can often really hurt performance; awk probably a better tool

Yes, both good points, example is still completed as per request. @terdon I specified line number as an issue, @Chris_Rands unlikely that fasta IDs would be formed to allow that kind of issue but point appreciated. FWIW, for manipulating fasta I would not use grep, awk, sed, there are plenty of good packages based around scripting languages, I would favour Bioperl. I am a lowly newcomer so couldn't comment this to OP.

Ah, so you did. Sorry, I missed the disclaimer and only saw the command. However, note that having having IDs be substrings of other IDs is absolutely common. For instance, the standard human genome fasta sequence has both chr1 and chr11.  That said, your answer has a much more serious issue: each line of rmfile.fa will be passed to grep as a pattern to search for, so any line containing that will be removed. This means that if sequence id3 is TAT, you will remove all lines matching TAT irrespective of their ID.

This is what samtools faidx is intended for.
It needs to be called twice.

First, to create a FASTA index (*.fai file):
samtools faidx input.fasta

Secondly, to filter the FASTA file with the help of the index:
samtools faidx -o output.fasta input.fasta ids…

ids… are the IDs to retain, as individual arguments separated by space. If you’ve got a file blocklist.txt with IDs you want to discard (one per line), you first need to invert this, after having created the index (using Bash syntax):1
remove_ids=($(awk '{print $1}' input.fasta.fai | grep -v -f blocklist.txt))

… and then we can invoke the command as
samtools faidx -o output.fasta input.fasta "${remove_ids[@]}"

What’s nice about samtools faidx is that it also works with BGZF-compressed FASTA data — i.e. fa.gz files (but only if they were compressed using bgzip, which ships with Samtools, and the index must be generated on the compressed file). However, it will always write FASTA. If you want gzipped FASTA output, you need to manually stream the result into bgzip:

samtools faidx input.fasta.gz "${remove_ids[@]}" | bgzip -c >output.fasta.gz

1 Careful, this syntax is generally unsafe for reading a command output into a variable; we can use it here only because we assume that we want to generate a space-separated list of identifiers; the generally safe way is a lot more convoluted.

Shortest way...
awk '!/Id2|Id3/' RS=">" ORS=">" fasta.fa

!= Not match the pattern
Note: Sequence doesn't need to be in a single line!

Firstly, you have to install biopython using the cmd pip install biopython.
Suppose you have two files id.txt and sequence.fasta. The solution file name is seq_rm.py which contains the following codes.
import argparse
from Bio import SeqIO
from Bio.SeqRecord import SeqRecord

if __name__ == '__main__':
    parser = argparse.ArgumentParser()
    parser.add_argument('--id_file', required=True,
                        help='file path of removable ids')
    parser.add_argument('--input', required=True,
                        help='path of fasta sequence')
    parser.add_argument('--output', required=True,
                        help='path of selected fasta sequence')

    args = parser.parse_args()

    with open(args.id_file, mode='r', encoding='utf-8') as file:
        ids = set(file.read().splitlines())

    records = [SeqRecord(record.seq, id=record.id, description='', name='') for record in SeqIO.parse(args.input, 'fasta')
               if record.id not in ids]
    SeqIO.write(records, args.output, "fasta")

Run cmd:
python3 seq_rm.py --id_file id.txt --input sequence.fasta --output selected.fasta

It will generate an output file named selected.fasta.

One way using awk without getline:
awk 'FNR==NR { a[$0]; next } /^>/ { f = !(substr($0, 2) in a) }f' ids.txt seqs.fa

I know you do not want a Python solution but nevertheless I would recommend it like others. Here is one way to do it using Biopython:
from Bio import SeqIO

ids_file = 'ids.txt'
input_f  = 'input_file.fasta'
output_f = 'result.fasta'

my_dictionary = {}
with open(ids_file, 'r') as f:
    for line in f:
        my_dictionary[line[:-1]]='value'

with open(input_f, 'r') as input_handle, open(output_f, 'w') as output_handle:
    for seq_record in SeqIO.parse(input_handle, "fasta"):
        if not seq_record.id in my_dictionary:
            SeqIO.write(seq_record, output_handle, "fasta")


text:
Are there certain alignment methods/tools which perform better with a high density of indels?
I have a set of experiments which should result in WGS reads with a high density of indels.
Question: Are there certain alignment tools/methods which perform respectively "better" with these reads? I'm not very interested in performance, but rather "better" in terms of dealing with the presence of many deletions/insertions.
Options I know of:

BWA  (I think there are special settings whereby I could run bwa
for this case, e.g. see this)
BBMap

try Stampy? ...

@Chris_Rands Have you seen this benchmarked against the others?

Which platform do you use? The link about bwa mem is about ONT platform. And what kind of indels are you talking about, small ones?

@b.nota Small indels, but a high density of these small indels

Minimap2 attempts to solve a lot of problems that are not properly dealt with by other mapping software. It's written by the same author as BWA, but designed to deal with situations where the reference and query sequences are quite different (e.g. for noisy long reads).

This is a great point, thanks! I've heard good things about this with PacBio reads

I will also keep this question open for a few more days; this is a great answer though!


text:
Which tools for differential expression analysis in scRNA-Seq?
I am starting to run analysis for differential expression in scRNA-Seq. Which tools are available for this kind of analysis?
Can tools for bulk RNA-Seq like DESeq be used for scRNA-Seq?

I recommend this page to learn about other tools

@Llopis thanks. I'd recommend also this curated repository for Single-Cell analysis: https://github.com/seandavi/awesome-single-cell

Yes, and there's also the online course, where it says that "Several computational analysis methods from bulk RNA-seq can be used"

According to Jaakkola et al. (ref. 1):

Based on our comparisons, DESeq and Limma without any modifications are not suitable for scRNA-seq data analysis, and yet, they have performed well in the context of bulk RNA-seq data.

They also provided a summary of the comparison between 5 tools (2 designed for Single-Cell, 2 designed for bulk, 1 general):

In my case I ended up using MAST (ref. 2, github).

Jaakkola, M. K., Seyednasrollah, F., Mehmood, A. & Elo, L. L. Comparison of methods to detect differentially expressed genes between single-cell populations. Briefings in Bioinformatics 18, 735–743 (2017).
Finak, G. et al. MAST: A flexible statistical framework for assessing transcriptional changes and characterizing heterogeneity in single-cell RNA sequencing data. Genome Biol. 16, 1183 (2015).

In a paper of Soneson & Robinson (2018) in Nature Methods (doi:10.1038/nmeth.4612), the authors concluded that "bulk" methods do not perform worse than methods developed specifically for scRNA-seq.

@b.nota Thank you! Can you edit my question to reflect also your comment?


text:
LRT or LRT-like test on cyclical (Sleep) data
I have RNA-Seq data from 4 time points (3 hours awake, 9 hours awake, 3 hours asleep, 9 hours asleep).  I'm interested in doing something similar to a LRT where genes are found to be significant if they changed at any point in the time course.
However, I'm not sure how to do this considering this is cyclical.  Running a LRT test in DESeq2 gives different results depending on which timepoint I choose as the beginning of the time course since I assume it doesn't calculate the last transition back to being awake (Or asleep depending on your reference).
Any ideas on an appropriate way to do this?
time = c(3,9,15,21)
ddsLRT <- DESeqDataSetFromTximport(txi.rsem.lrt, sampleTable, ~time)
ddsLRT<- DESeq(ddsLRT, test="LRT", reduced = ~1)

I think the problem is that you have time as a linear value rather than a factor. While this naively makes sense (after all, you have to go through 3 hours awake before you can get to 9 hours awake), what you end up doing is fitting a single coefficient to each gene of "linearly changes with time". What you want instead is factor, so you can have a system where things like per and clk can cycle:
time = factor(c(3,9,15,21))

The rest is then correct. Note that LRTs are not asking, "is there a change at any time point?", since there will likely be DE genes according to this that are not DE at any single time point. Rather an LRT such as this is asking, "What are genes that vary significantly over time?" There's a subtle different between those two questions. The former question is actually answered by a bunch of pairwise comparisons (taking the union of the DE genes).

Thanks for this!  I actually tried this last night after posting and it seemed to work, but I wasn't 100% sure it meant what I thought it meant until your post.  Much appreciated.


text:
Filtering bam file with bedtools intersect using strand specificity
I have a small sample of records from a bam file. I obtained them by intersecting with a bed file not enforcing strand-specificity. So I get them all back upon performing the same intersection again:
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 | samtools view
MN00481:17:000H2FH5N:1:11104:12480:15984    16  III 4499483 42  28M *   0   0   CAGGAACACGCACATCGTCGATTTTCAG    FFFFFAFFFFFFFFFFFFFFFFFFFFFF    AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:28 YT:Z:UU
MN00481:17:000H2FH5N:1:11103:9925:11522 16  IV  2388448 42  30M *   0   0   CACTGGAACTGGTCCGATCTTCAAGCCGGT  FFFFFFFFFFA/FFFFAFFFFFF/FFAFAF  AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:30 YT:Z:UU
MN00481:17:000H2FH5N:1:22109:11939:5243 16  IV  4749336 42  29M *   0   0   ATGTGGATGCGGTCAGCTTGGGCACGATG   FFFFFFFFFFFF/FFFF6FFFFFFFF/FF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:22104:17385:19431    16  V   4443064 42  29M *   0   0   CAAAATGTGAACAGGAACGTTTTCATGAG   F/FFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:23105:22058:12968    0   X   981524  42  29M *   0   0   CCATCGATTCTCATTTTAATGCTGGTGCT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU
MN00481:17:000H2FH5N:1:12104:7867:7500  16  X   4214381 42  29M *   0   0   GCTGCAGGAATGACGGCTCTAGCAAACGT   FFFFFFFFFFFFFFFFFFFFFFFFFFFFF   AS:i:0  XN:i:0  XM:i:0  XO:i:0  XG:i:0  NM:i:0  MD:Z:29 YT:Z:UU

Now I would like to enforce strand-specificity, but in either case, none of the records pass the filter:
Same strand (-s):
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -s | samtools view
$ # nothing here

Opposite strand (-S):
$ bedtools intersect -a sample.bam -b features.bed -wa -u -f 1.0 -S | samtools view
$ # nothing here

What is the correct way of doing this?
Edit
Bedtools version:
$ bedtools --version
bedtools v2.27.1

My bed file has strand information:
$ tail features.bed
X   17617826    17619458    "WBGene00015867"    +   .
X   17622777    17640743    "WBGene00019211"    -   .
X   17650177    17651594    "WBGene00044771"    +   .
X   17662457    17665357    "WBGene00003744"    +   .
X   17671028    17680790    "WBGene00020732"    +   .
X   17680910    17682294    "WBGene00044772"    -   .
X   17695231    17701744    "WBGene00017625"    +   .
X   17710485    17711674    "WBGene00019189"    -   .
X   17714455    17714832    "WBGene00007069"    +   .
X   17714967    17718717    "WBGene00007068"    +   .

Your bed file switch the score and strand columns. Strand should be column 6, score column 5. For example take this bed file with the 5th and 6th columns switched - zero results for strand specific intersection b/c strand is actually '.'
$ cat wrong.bed
chr1    30000000        31000000        noname  +       .
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 \
    | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b wrong.bed -wa -u -f 1.0 -s \
   | samtools view | wc -l
0

After fixing the bed file:
$ cat right.bed
chr1    30000000        31000000        noname  .       +
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 \
      | samtools view | wc -l
166
$ bedtools intersect -a mm10_500k_pe.bam -b right.bed -wa -u -f 1.0 -s \
      | samtools view | wc -l
83

Your score and strand columns are swapped. If you correct that then you'll get stranded results.


text:
plotMA function issue: can not circle genes
I am trying to make MA-plot for my bulk-RNA-seq dataset and experiencing issues with it. I just copy-pasted the code that can be found here:
https://f1000research.com/articles/4-1070/v2
which is:
library("DESeq2")
plotMA(resLFC1, ylim=c(-5,5))
topGene <- rownames(resLFC1)[which.min(resLFC1$padj)]
with(resLFC1[topGene, ], {
  points(baseMean, log2FoldChange, col="dodgerblue", cex=2, lwd=2)
  text(baseMean, log2FoldChange, topGene, pos=2, col="dodgerblue")
})

However, I am not getting the right plot with the top gene circled, but the following one:

Maybe it is a bug in the plotMA function, not sure. How can I reproduce the figure from the publication?

Are you sure that resLFC1 is created as stated in the article ? There seems to be a lot more points than in the article. Also, do you get any error or message?

Are you sure the weird double quotes (“ versus ") are not the problem?

I corrected them, of course. Rstudio complains about it right away.

Your plot shows all genes as DE (red points). That’s suspicious.

It's likely that the fold-change of the topGene is greater than 5. Consequently, the circle is being drawn and the text produced, but they're outside the bounds of the plot (the bounds don't update if you add a circle and text label). Look at the log2FoldChange column of resLFC1[topGene, ] and increase the bounds to encompass that.

license: Creative Commons - Attribution Share-Alike - https://creativecommons.org/licenses/by-sa/4.0/


text:
Highly variable genes analysis by DESeq2
In the following article. There are two ways of getting highly variable genes:
resSig <- subset(res, padj < 0.1)
head(resSig[ order(resSig$log2FoldChange, decreasing=TRUE), ])

Getting genes that are statistically significant, then ordering them by the log2FoldChange which is

The column log2FoldChange is the effect size estimate. It tells us how much the gene’s expression seems to have changed due to treatment with dexamethasone in comparison to untreated samples.

And the second one:
topVarGenes <- head(order(rowVars(assay(rld)),decreasing=TRUE),20)

I am struggling to understand what are the differences between the two ways of getting the most variable genes. What is the difference?

The first method (resSig$log2FoldChange) is for getting the most differentially expressed genes. These are the genes with the biggest differences between specific groups you pre-define.
The second method (rowVars(assay(rld))) is for getting the highly variable genes. These are the genes most variable across all samples regardless of which samples they are.


text:
Why do NEB adapters have non-complementary sequence?
I am making some infographics of different library prep types and noticed something weird about NEB's Ultra II adapters.
Molecule 1 is the desired product of the adapter ligation step. * indicates a phosphorothioate bond and - indicates a normal bond. | indicates complementary sequence through hydrogen bonding.
Molecule 2 is the product after treatment with the USER enzyme.
Molecule 3 is the result of melting molecule 2 and hybridizing the NEBNext Universal PCR Primer for Illumina.
My question is: Why does Molecule 3 have a large non-complementary region on the 3' end of the original molecule? Is that just to force Molecule 1 to have an open loop conformation so that the USER enzyme can access it better? Or maybe it is something to do with PCR kinetics?
I understand that PCR synthesizes in the 5'->3' direction, so Molecule 3 only forms with the original adapter-ligated molecules before PCR.

Module 3 does not exist during PCR amplification. Here is how NEBNext adaptor and multiplex PCR primers work.
During the 1st PCR cycle, only the P7 primers are used to bind and amplify ligated DNA, generating the complementary strands that contain the binding sites for P5 primers (Universal in this case). During the second cycle of amplification, the P5  primer hybridizes to the DNA strand generated during the first cycle and start to amplify. Please refer to the following figure for more details.

The original point of the non-complementary sequence is to allow formation of a Y-adapter after USER treatment, since that's what's needed for the Illumina chemistry to work. In reality the bottom strand may not be there in step 3, since USER has already digested one strand (creating the hairpin in step 2) and size selection has been performed. Regardless, it'll get removed once cluster formation is performed.


text:
How long does it typically take to call peaks via `phantompeakqualtools`?
I'm running the script run_spp_nodups.R from phantompeakqualtools. The input is a pair of BED 3+3 ("tagalign") files that are roughly 300M, and I'm using 8 threads. How long would you expect it to take? I am happy to provide other info to the extent that I can manage -- post me a comment!

Is there a reason you want to use phantompeakqualtools? It's not exactly commonly used.

If I were building a pipeline from scratch, I probably wouldn't. I am using it as part of Anshul Kundaje's open-source pipeline.
https://github.com/kundajelab/chipseq_pipeline
I'd rather not go in and start rewiring things on my own if it's running properly -- but it's taking long enough that I suspect that something is off. Hence the question.

Is there a specific reason that people avoid it?

It's a weird hodge-podge of cryptically written and poorly documented scripts that may or may not give reasonable output. It's basically an artifact of a time when ENCODE was producing a LOT of data and no one had any clue how to reasonably QC all of it. So the NSC/RSC metrics (the main point to phantompeakquals) were created. In hindsight, they're not terribly informative, in part due to being incredibly sensitive to settings. Just do some basic QC and call macs2. The pipeline will do just that anyway.


text:
How can I use UpSet plot for expression/RNA-seq data
I have a data.frame of 27 000 x 17 and I want to feed it into UpSetR program, but getting the error:

Error in start_col:end_col : argument of length 0

My data.frame looks like that:

And the upSetR code that I am using:
upset(df, nsets = 4, number.angles = 30, point.size = 3.5, line.size = 2,
mainbar.y.label = "Gene Intersections", sets.x.label = "Genes Per Group",
text.scale = c(1.3, 1.3, 1, 1, 2, 0.75))

What am I doing wrong here? Any suggestions would be greatly appreciated.

Update

I found here a similar issue:
https://plus.google.com/109635504520741246387/posts/Etujrp6sU1Q
They say that the input should be binary. But then, how can I calculate this binary matrix from the rlog-data? I am not sure...

UpSet plots are in many ways similar to Venn diagrams, in that the represent some sort of metric for group overlap. What I expect you want to do is get the number of DE genes in various comparisons and put those numbers into an UpSet plot.

upsetR takes as input a data.table in the long format. I found this example/code very helpful when I first started using upsetR.

I need to get the binary matrix to be able to feed it into the UpSetR. Right now I have a matrix where rows are genes and columns are groups. Somehow for each gene I need to set a threshold value to convert it to binary ‘expressed/not expressed’ probably. Groups are vectors of numbers, so I need to find their intersection somehow. I do not know which software or R package do that: I would not want to do it manually.

I am pretty confused but the problem is rather definite: I need to represent groups, their similarity, as an UpSet plot.

I think, it would be good to rephrase your question to match your ultimate goal something like "How can I use UpSet plot for expression/RNA-seq?". Say what is you sample and why you want to do the UpSet plot.

You need to generate a list of gene ID vectors to feed into upset. In the example below, I am selecting observations (genes) which have an level above the median level. You could also choose only genes above the 3rd quartile, etc. You then need to convert these from a list of character vectors to a presence/absence matrix. UpSetR has a fromList function that does just this.
In future, be more specific about what you want to achieve and supply example data for people to work with.
library("UpSetR")

tmp <- mtcars

tmp$car <- rownames(mtcars)
list <- lapply(colnames(mtcars), function(x) {
  tmp[tmp[[x]] > median(tmp[[x]]), "car"]
})
names(list) <- colnames(mtcars)

upset(fromList(list))

Created on 2019-10-17 by the reprex package (v0.3.0)

I wrote a script to create the binary data frame with my differential expression data. If it is of any help to anyone:


text:
Performing t.test to obtain p.values for each gene
I have FPKM values for Normal(80 cases), Group1(135 cases), Group2(147 cases), and Group3(102 cases) for 160 genes (Groups were obtained by clustering tumor data).
I want to perform t-test (two-sided, unpaired, unequal variance) between Normal and Group1, Normal and Group2, and Normal and Group3 such that I get p.value for each gene.
The t.test function in R returns single p.value and not gene-wise (I don't want to do t-test in Excel).
Which package/function should I use to do this in R?

You can simply apply() t.test() to your matrix. In general, though, I expect lmFit() from the limma package will be simpler to use.
Here is an example with t.test():
m = matrix(runif(100, max=c(rep(1, 50), rep(1.5, 50))), nrow=10) # make some data
apply(m, 1, function(x) t.test(x[1:5], x[6:10])$p.value)

Neither t.test() nor lmFit() should be expected to produce the same results as excel, since excel is inferior in every way and should never be relied upon for any statistical analysis (plus lmFit is doing something much much more complicated than a simple t-test).
As an aside, FPKMs are highly problematic for use in statistics. If possible, I strongly encourage you to get either per-gene counts or per-gene estimated counts and to use those in limma or another differential expression package instead, as the results will be much more reliable. I presume, though, that you're using TCGA data, so you may be stuck with FPKMs.

I totally agree with Devon's aside comment. FPKM values are not suitable for (sound) statistics and only informative for heatmaps etc. For proper statistics you should use count data, all statisticians agree on that.

@b.nota Source?

@SmallChess, for example paper about edgeR, but also the more difficult stuff of Gordon Smyth and read the discussions on bioconductor support.

I tried t.test()Group1.ttest <- t.test(data.Normal, data.Group1, "two.sided", paired = FALSE, var.equal = FALSE)but the function does not give me p.values row-wise. The result that I obtained for above code> Group1.ttest$p.value [1] 5.632393e-48I tried lmFit() fac <- factor(c(rep(0,80),rep(1,135))) fit = lmFit(Normal.Group1, design=model.matrix(~ fac))   e = eBayes(fit) tab<-topTable(e, number=nrow(Normal.Group1))  but the p-values obtained are different from the ones obtained when I do the t.test on the same data in Excel. I tried log2(Normal.Group1) but still the results differ.

The results should differ, excel should never be relied on for anything important. lmFit is also a more complicated computation than a simple t-test. Please see my updated reply.

Is a  parametric t-test really suitable, either?

Depending on the input, yes. Limma is basically a fancy t-test.


text:
What's "signal" in gene expression?
I'm exploring gene expression on public available databases. RefExa is one of the databases I've been studying.
When I searched something there, the website reported the expression value unit as "signal". See screenshot where I highlighted:

I'm not able to find anywhere what that really mean. TPM? RPKM? Raw counts? Normalizated for gene length?
My question is that how to interpret gene expression "signal"? What does signal 100 mean?

Looking at this resource, these data are microarray data rather than RNAseq data. Thus the "signal" in question is the intensity of florescent light being emitted from the the probe associated with your gene.  The relationship between intensity and quantity of RNA is different for every probe because each probe will have a different affinity for complementary fragments dependent on the sequence in question (not least the melting temperature).
Affy intensities are usually presented as log values going up to around 12, but the raw number is next to meaningless and usually only means something as a normalized (e.g. RMA) log fold change.
I've never come across RefExA before, and it appears to be undocumented and unmaintained. I wouldn't rely on it for usable data if I were you. A more modern database of expression data is the EBI Expression Atlas. It presents data from a wide range of different studies, processed and normalised on a unified pipeline and expressed as TPMs.


text:
BioMart returns one more item
I am supplying 30 ensembl ids into BioMart, but somehow it returns 31 result. What is going wrong here?
Here you can see the code:
library('biomaRt')
mart <- useDataset("mmusculus_gene_ensembl", useMart("ensembl"))
genes <- rownames(mat)
G_list <- getBM(filters= "ensembl_gene_id", attributes=
c("ensembl_gene_id","entrezgene",
"mgi_symbol","description"),values=genes,mart= mart)

genes list is the following one:
[1] "ENSMUSG00000097971" "ENSMUSG00000097582" "ENSMUSG00000061808" "ENSMUSG00000097206"
[5] "ENSMUSG00000022425" "ENSMUSG00000094708" "ENSMUSG00000095937" "ENSMUSG00000032484"
[9] "ENSMUSG00000060550" "ENSMUSG00000073409" "ENSMUSG00000026822" "ENSMUSG00000056071"
[13] "ENSMUSG00000035299" "ENSMUSG00000045932" "ENSMUSG00000006154" "ENSMUSG00000032496"
[17] "ENSMUSG00000040809" "ENSMUSG00000032690" "ENSMUSG00000055322" "ENSMUSG00000066443"
[21] "ENSMUSG00000005089" "ENSMUSG00000073489" "ENSMUSG00000005800" "ENSMUSG00000034855"
[25] "ENSMUSG00000029417" "ENSMUSG00000029561" "ENSMUSG00000021843" "ENSMUSG00000020108"
[29] "ENSMUSG00000037849" "ENSMUSG00000066861"

As you can see I have here 30 genes, but BioMart returns:
> G_list$mgi_symbol
[1] "Slc1a2"  "Mmp8"    "Eps8l1"  "Ddit4"   "Ktn1"    "Enpp2"   "Lcn2"    "Cxcl9"
[9] "Oasl2"   "Ngp"     "Ltf"     "Oas2"    "Cxcl10"  "Mid1"    "Ifi206"  "Chil3"
[17] "Ifit2"   "Tns1"    "S100a9"  "H2-Q7"   "H2-Q7"   "Ttr"     "Gm10163" "Oas1g"
[25] "H2-Q6"   "Ifi204"  "Gm10359" "Gm12671" "Gm26699" "Gm26527" "Gm26917"

If you print G_list, you'll see this half way down:
1 ENSMUSG00000060550      15018      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]
2 ENSMUSG00000060550     110558      H2-Q7      Mus musculus histocompatibility 2, Q region locus 9 (H2-Q9), mRNA. [Source:RefSeq mRNA;Acc:NM_001201460]

So there are two entrez ids that map to the same Ensembl id. One of them maps to an alt locus.
If you exclude "entrezgene" from the query, you'll get 30 results.