text:
RefSeq transcript naming - transcript vs gene
I have a dataset composed on numerous RefSeq transcripts with names of the form NM_###### - essentially standard annotated mRNA accession numbers.  I notice that numerous transcripts which have different sequences (different start and stop sites in a gene) all have the same NM_#### accession number.  The interpretation is that they are all from the same gene, but this makes no sense to me.  Why should different transcripts with different sequences have the same name?  I understand this identifies their gene of origin, but they are still different transcripts.
Here is an example from the .bed file I'm using.  Note they have different start and stop locations and even have different strand notations, but all have the same name.
chr     start        stop       name            score   strand
chrX    3441730     3443690     NM_001085545    0       -
chrX    3700232     3702192     NM_001085545    0       +
chrX    3750925     3752885     NM_001085545    0       -
chrX    31336368    31338328    NM_001085545    0       +

Can you provide a few examples? If multiple entries have  same NM_### accession but different sequences, could that be because they are different versions? You can view the entire version history of any NM_### accession on the NCBI pages.

What you're describing is almost certainly a mistake in your data. Different sequences don't have the same refseq IDs. Please give us an example.

I added an example of what I see in the .bed file.  Note that I misspoke a bit - I don't have the actual nucleotide sequences, but these all have different genomic start and stop sites, so they are effectively different sequences.

Are they the exons?

The current ensembl entry for that transcript is the only entry for that gene, it has a single exon spanning 3.70 to 3.72 Mb.  So what is that entry at 31 Mb supposed to represent?

Searching for the transcript name in ENSEMBLE, I see the location and strand are the same as the second line in my example (almost the same:  start = 3700233 - one different than mine, stop = 3702192, forward strand).  It is also noted as having only one exon.  I don't understand where the other 3 entries are coming from.

We can't tell you that. Only the person who made the bed file can. What is the bed file supposed to be showing? As I show in my answer below, the sequence NM_001085545 can be mapped extremely well to multiple regions of the mouse genome. So why many of those matches are included in the bed file is a question for the person who made the file.

No, this is the same sequence. A single RefSeq accession will correspond to one sequence. However, nobody said that sequence needs to be unique in the genome. The specific sequence you are talking about, NM_001085545.1, is a mouse transcript which can be aligned to many different regions of the mouse's chromosome X.
This is your sequence:
>NM_001085545.1
agagaagtctaggcaggactggcagactacagcaatgcctcaaccaagaagcctctgcga
gagcttgttctaggaattttcatcaatggggcttttagtcagcagggtcttgagatgcag
ggattccagtctgttagagccacagccagaagccatagccggagccagctacattcctgg
cagtcgcaagcgaaaaagaaacagtttggaggagttggcaacaagttctaatgttcatgg
acctcaaaaccagggaatgtatccacatcaagttctcaactacatctactggaaaagggt
taagatctcatctaatgatgcttatcaaaacttatttttggacgggcatgatagcgacat
taaaatccgtgctctgggaaaaacatggtgtttacacaaagtatttttatgtcagtcagg
ctactttgctaacattctcaaaggtacttggagagaatcacaccatggtgttataaatct
gatcattaagaatgaggatattgatacccgatctctgcattttgtgtttggtgctttgta
cacggatgcggatttgtcaataacacctctggaagttcctcaagttttggcagcagcatg
cctgcttcgggtggatcaagtaattcagcagtgtgaaggaatcatgaaagaaactatcaa
caggaacactgtgtgctcctattatttggcagcagaaacctatagattaaaagctgtaaa
gacgagatgctttgaatggcttctttgcaatttgatggtacatccaagtgtggcacttta
caaggaagtagacttgaagttgatgtatcttctagcactgtcttctgacttactagtcat
gcaaaaggagattgatgtatataccacactaaaaatatggatgttcctttatcttaatcc
atgctggaacggaaccatgaaacagcttttacaacacgcaaacaactggctttccaccca
catggcatatgttgataacatcagttttcttgaaagtgaagaaggactaatatttcaacc
agtgtttaaaaagctgagatttcagcacattatctgtgacttgacttccacaactattct
tgaacaagatcgactaatacctatggcatggttgtcacccatttacaaacaacagtggtt
gactttgctgcgaacacaagaatatggggtaattggaccacaagttatcaatgaacaaga
acttgaagaatgcaccatgaggtgtggtacaatgatccccaaggatggaagatatacttg
gaagtggtcagttggacgacttggctttcctttacgtgtgacctttaccaggcagtgtgt
aattttaaggcaacggtgtcagaggtgtgatggttctgcttgccacaaccatatccgaaa
tgtcatattcagaataactttggtgtgttttgattccaacaaaagagtaactttcagaaa
gacaacaggttataaaatcctcacctttgaatataaggaggagcaaattgtaatgaaatt
ggatagtgatgttctaaccttccctatgtgtatattctgcaatttcctttttgtaaacct
aggaaatgcagaaaacaagtaatcttatcaattattagactgtaagtattttacaagtat
tgaatggttcatttaatatgaaggctgtgtaatgactacagtgaagaagaccaaaaacaa
attcaatattctttggagtctagtaccaccagaagtgactgcttccctccatcaccttat
ttccatgtatttctgaagaaaacctttacttccaacaaagtatttgtcaaatcaaataaa
cctaatcttctgatttaaatttaccatcgagatccaagcttggaatggaaataagagact
atatttgcaaggtaaaacatctcaatcaatttgttaagtagatttggccctagctaaaat
tgagacaaaataaaatcatgccattgttgaaactgttgtg

If you BLAT it against the mouse genome using the UCSC BLAT tool, you will see many hits, most of them almost perfect matches and some of them overlapping the regions given in your bed file:

Your bed file matches the following hits in the image above:
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3441731   3443690   1960
NM_001085545.1  1960     1  1960  1960   100.0%  chrX   +     3700233   3702192   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   -     3750926   3752885   1960
NM_001085545.1  1954     1  1960  1960    99.9%  chrX   +    31336369  31338328   1960

I don't know why your bed file has all of these various matches for NM_001085545.1, since you haven't told us what the bed file is supposed to be doing. A bed file is just a collection of regions. For all I know, that bed file is just "a few hits for sequence NM_001085545.1 in the genome". There is no reason to assume that each region in a bed file corresponds to a gene.
So there's nothing wrong with your file, you're just dealing with a non-unique sequence from the mouse genome. And you don't have multiple sequences with the same accession, you have one sequence, one accession, but multiple hits for that sequence in the genome.


text:
Fast processing of fastq data
I am trying to write python script for customized filtering for fastq file (size >3 GB). My proposed script is as follows:
def filtering(read):
    time.sleep(0.1)

    if len(read) >= 15 and \
       len(read) <= 30 and \
       np.mean(read.qual) >= 30 and \
       sum(i < 20 for i in read.qual) <=2:
           flag1 = True
           return  flag1

pool = multiprocessing.Pool(5)
pool.ncpus = 4
for file in os.listdir('/moRNA_data/trimmed_fastq'):
    if file.endswith('.fq'):
        fastq_file = HTSeq.FastqReader( os.path.join('/moRNA_data/trimmed_fastq',file))
        name = str(file.split('.')[0]) + '.fastq'
        file_out = open(os.path.join('/moRNA_data/filtered_fastq',name), "w" )
        for read in tqdm(fastq_file):
            flag1 = pool.apply_async(filtering,read)
            if flag1:
                read.write_to_fastq_file( file_out )
pool.close()
pool.join()

The problem with this script is that it works well with small file like upto 100-200 MB but when it goes to directory having multiple sample each having size more than 3
GB, it works well for first file but after that it hangs so severely that I had to restart the server [ I can't even kill the process]. I need suggestions to improve this script to make fastq file processing faster. Platform details is as follows:
OS : Ubuntu-18.02
Architecture:        x86_64
CPU op-mode(s):      32-bit, 64-bit
Byte Order:          Little Endian
CPU(s):              12
On-line CPU(s) list: 0-11
Thread(s) per core:  2
Core(s) per socket:  6
Socket(s):           1
NUMA node(s):        1
Vendor ID:           GenuineIntel
CPU family:          6
Model:               158
Model name:          Intel(R) Core(TM) i7-8700
CPU @ 3.20GHz
Stepping:            10
CPU MHz:             800.560
CPU max MHz:         4600.0000
CPU min MHz:         800.0000
BogoMIPS:            6384.00
Virtualization:      VT-x
L1d cache:           32K
L1i cache:           32K
L2 cache:            256K
L3 cache:            12288K
NUMA node0 CPU(s):   0-11

Thanks.

Make sure that filtering() actually returns something, though that's likely not the root of the problem.

@DevonRyan just to clarify, tqdm only reports time in such cases it does not do a list conversion as you indicated- you can easily test with with an infinite generator https://pastebin.com/NJvwf5DT

Take out the time.sleep() call for starters. and I think you're missing an else: flag1 = False clause. If you still have issues, remove the multiprocessing code to assist your debugging

Did you try it on just one CPU or two just to test if the problem is with the threading and other applications on the server? Also how much RAM and memory is available?

@llrs : This codes works amazingly with small files. As I have mentioned in questions. I also tries with Pool(2) but same server hangs after processing of one file very severely.

@Chris_Rands : My aim to ignore those reads which do not satisfy the given conditions. So I don't think it change something even if I write else statement. For tqdm approach, my only aim is to check how fast this script is checking sequences. When I run this script then for the first file, I saw that it is processing 8000-10000 reads per second (very fast) and process a fastq file of size 3 GB in few seconds. But after that it's block everything. I can't even kill the process. Thanks.

@DevonRyan: I checked that too by printing it on small test file. But it prints only multiprocessing object details rather than its actual value. I also tried flag1.get but gets object description again. Again On small file, this script is OK. Many thanks for positive EDIT in question to make it more informative.

you can ignore my comment about setting flag1 = False, it was an issue that the previous edit caused changing the indentation. Now (after my edit) the function returns None (which is falsly anyway) without an else-clause. you should still remove the sleep() call but i think your issue is due to the multiprocessing, strip that part out and see if it works

@Chris_Rands : Sorry if I am not able to understand your suggestions. Multiprocessing is the main part of program and my main goal is to split filtration of millions of reads in multiple process to make it faster. If I remove it then it will be a simple for loop script and definitely will work but extremely slowly. So what will you recommend in place of multiprocessing.

it's unclear to me if you want to multiprocess per file or also within each file (per read). the per file is multiprocessing is much simpler. personally, i'm not a big fan of python's multiprocessing module so i would solve this (per file) via a pipelining wrapper like snakemake or ruffus

@Chris_Rands: It is not multiprocessing per file but it is within file. If you check the argument I am passing in pool.apply_async(filtering,read), I am passing a read sequence to a function. Read sequence I am getting using for loop. I may have little knowledge about per file application of multiprocessing. As far as I know, I can use it on very heavy file to make filtration task faster to split in many parallel processes.

If you can't kill the process and your only option is to restart your server, it seems that you exceed your available RAM. So you could monitor the RAM usage at first and check how to clear the memory.

@Lot_to_learn As detailed in the multiprocessing API docs, Pool.apply_async returns an AsyncResult object, not the result of the call directly, so I don't think your code is doing what you expect.

If this is not just an exercise for you, I would suggest to use already existing tools. bdduk.sh from bbtools should have all the options you need. And btw: You cannot use np.mean to calculate the average quality. Because the phred scores behave logarithmic  and not linear.

Note that str(file.split('.')[0]) is likely not doing what you want when a file name contains more than one dot. And it seems to me that str is unnecessary: you already have strings. You might be interested in os.path.splitext.

@bli: That's I have taken care of. All files names are Sample_ID.fastq . So In all the files there is only one dot. For str you are right. I can remove that.

@mbrig : I also realize the same that applying if condition to Pool.apply_async object may be the reason to hang system because this may set an infinite loop.

@finswimmer : Thanks for suggesting bbduk. I have check. I have found quality trimming based on mean quality and based on length. There are still two consition for which I have to write script. Filter reads having more than 2 bases under quality score 20 and filter reads with unique sequence having counts less than 10 in fastq file.

Knowing a high-performance language will make a huge difference here. See the example below in C. I haven't tested, but it should be easy to modify for your purpose.
C++, Rust, Go, Nim and Julia can be as fast.
// Download https://raw.githubusercontent.com/lh3/minimap2/master/kseq.h
// Compile with: gcc -O2 -o myprog this-file.c -lz
#include <zlib.h>
#include <stdio.h>
#include "kseq.h"
// STEP 1: declare the type of file handler and the read() function
KSEQ_INIT(gzFile, gzread)

int main(int argc, char *argv[])
{
    gzFile fp;
    kseq_t *seq;
    if (argc == 1) {
        fprintf(stderr, "Usage: %s <in.seq>\n", argv[0]);
        return 1;
    }
    fp = gzopen(argv[1], "r"); // STEP 2: open the file handler
    seq = kseq_init(fp); // STEP 3: initialize seq
    while (kseq_read(seq) >= 0) { // STEP 4: read sequence
        int i, sum, low;
        if (seq->qual.l == 0) continue; // no quality
        if (seq->seq.l < 15 || seq->seq.l > 30) continue;
        for (i = sum = low = 0; i < seq->seq.l; ++i) {
            int qual = seq->qual.s[i] - 33; // assuming Sanger quality
            sum += qual;
            if (qual < 20) ++low;
        }
        if (sum < seq->seq.l * 30.0 || low > 2) continue;
        printf("@%s\n", seq->name.s); // output can be made faster
        puts(seq->seq.s);
        printf("+\n");
        puts(seq->qual.s);
    }
    kseq_destroy(seq); // STEP 5: destroy seq
    gzclose(fp); // STEP 6: close the file handler
    return 0;
}

To apply the program to all *.fq in a directory (again, not tested):
mkdir -p filtered
ls *.fq | xargs -i echo ./myprog {} \> filtered/{} | sh

You can also replace the final sh with parallel, but the above should be fast enough.

:Thanks for your response. Till now I have done all my analysis in either python or R. I have never touched C/C++. Although I studies those language but I left using this before 8 or more years. But still I think I should try your way [although it is very much time taking for me to start heavy text file processing in C/C++]. For time being, Do you have any python alternative solution for this problem?

@Lot_to_learn The inner loop mimics your script. It won't take you much time to modify it for other purposes even if you don't know C well. For tasks like this, C will be times faster than python whatever you do.

source: Stack Exchange

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

sort: votes


text:
How drastic would the result be if I use fasta or reference assembly from ucsc and gtf from gencode?
There are difference annotation file for UCSC and gencode.
But if I use the reference assembly from UCSC and the GTF from Genocode or vice versa would my downstream results would be wrong?

Do you want to include older UCSC/gencode releases or only the most recent ones? They should now be fully compatible.

so i used gencode.v21 as my gtf file and not sure about the reference assembly,where i have used hg38 but im not sure its from gencode or ucsc .Is there a way to find out the what reference assembly im using from the reference assembly ?

Never use genomes or annotations from UCSC, they're poorly versioned and only recently with mouse and human have they even included all of the contigs. For fasta/GTF files from early in the GRCh38 release, you can tell whether you're using UCSC or Gencode by the presence/absence of _random contigs, which will only exist for UCSC. These were mostly later split into the actual contigs, so recent download from UCSC should more closely match what you find at Gencode/Ensembl. Further, that time predated UCSC beginning to adopt Gencode's vastly superior annotations, so if your GTF file has instances where the same gene ID is on either multiple strands or multiple chromosomes (this is obviously biologically impossible) then you have a UCSC GTF file.
In general, with early GRCh38 releases your only real issues will be with _random contigs, which are a minority of the genome and don't have all that many genes. But really, you should be keeping track of the sources of your files and ensuring that they're compatible.
Update: I should expand a bit on my "Never use genomes or annotations from UCSC" comment. In point of fact the genomes themselves aren't so terrible. Early on UCSC had the bad habit of concatenating contigs together into _random "chromosomes", but they seem to have mostly kicked that habit as of late. Note, however, that there are no versions for their genomes. Since reference genomes continue to get updates over time (mostly through the addition of patches) the lack of actual versions means you have to manually check if a recently downloaded file matches what may have been downloaded either previously or by someone else. This has obvious consequences when it comes to reproducibility. The same issue occurs for annotations from UCSC, but they have the additional problem of historically having biologically incoherent concepts of genes. That is, they will contain the same gene in multiple places with multiple orientations, which will break many tools in both obvious and completely unclear ways. For example, DEXSeq will simply break with an error message if given a UCSC annotation, since they break biological plausibility. If you were to use these annotations files with deepTools, you wouldn't get an error message, but the resulting output would be only partial, due to the biologically impossible annotation effectively corrupting most obvious ways of storing annotation data in a data structure (i.e., you can no longer treat IDs as unique). This could have downstream ramifications on biological interpretation of results.

to make sure if im have gencode so i used this "cat gencode.v21.annotation.gtf | grep "_random" no instance of that word .

You could also just head the file. Gencode files tend to start with a few comment lines (granted, given what the file is named, the odds of it being from Gencode were very high to begin with).

##description: evidence-based annotation of the human genome (GRCh38), version 21 (Ensembl 77)
##provider: GENCODE
##contact: gencode@sanger.ac.uk
##format: gtf
##date: 2014-09-29   im sure now im in right path

@krushnachChandra not really relevant, but just so you know, grep can read files. You don't need cat file | grep foo, you can always do grep foo file.

sure i was just trying to see if i haven;t messed up my annotation file .I will surely keep that in my mind. "grep"

Devon, could you elaborate on your advice never to use UCSC genomes. Never seems a bit extreme. Can you explain what they're missing in a bit more detail?Is there any reason not to use them for human reference? Yes, they may have *random in them, but is that such a problem?

*random is the least of the problems. I'll elaborate more in my answer.

Thanks, that makes sense. I have used the (human) UCSC genome sequences extensively and haven't had any issue with them although I admit that the lack of versioning is a real pain. I was afraid you were suggesting the actual sequences were somehow of lesser quality.

The main difference is in the way the chromosomes are named - UCSC uses the "chr" prefix (so chromosome 1 is "chr1") while in gencode the "chr" isn't used (so chromosome 1 is just "1"). Depending on your use case, this can obviously cause problems - if you're trying to match a locus (e.g. from gencode 1:1000002) between them, whatever tool you use is going to be looking in your aligned data for "1:1000002", but in your aligned data, it'll be named "chr1:1000002", so it won't match the two up.

Both use the chr prefix for the most recent human and mouse releases.

so its with the annotation file my results would vary  but not much with the references?

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


text:
How to reverse complement the DNA sequences for given inverse/reverse coordinates?
I have the series of coordinates in id.txt file, whose coordinates sequences are in genome.fasta file.
The coordinates of id.txt file is shown below,
Contig3:15-7
Contig2:5-10
Contig1:12-3

The genome.fasta file is shown below,
>Contig1
AAGGCCATCAAGGACGTGGATGAGGTCGTCAAG
>Contig2
ACGAAGCGCTCGCCAAGGCCGAAGAAGAAGGCC
>Contig3
GCTGCGGCGCTGATCCTGGCGGCCCGCGCCGAG

I have used the following codes to extract the sequences from genome.fasta based on the coordinates in id.txt by using the following command xargs samtools faidx genome.fasta <id.txt > result.fasta
It has extracted only Contig2:5-10sequence. Because, the coordinates of Contig2:5-10is in proper order. But Contig3:15-7 and Contig1:12-3 coordinates are in inverse/reverse order, So samtools could not fetch those sequences.
I need to extract those sequences from Contig1:12-3 and Contig3:15-7 coordinates and also I want to reverse complement them. I have  larger fasta file and plenty of coordinates sequences to be extracted, So I need to automate this process. So, please help me to automate this process.
The expected output is
>Contig3
GATCAGCGC
>Contig2
AGCGCT
>Contig1
CTTGATGGCC

Thank you in advance.

Welcome to the site. Do you have any preference about the tool to be used? What else have you tried to set the coordinates in the right order? Usually this kind of problems are solved by working on every tiny problem one after the other.

I have used linux command the following command to keep the coordinates in right order,                                                                                "awk '$2 > $3 {printf( "%s\t%s\t%s\n", $1, $3, $2); next} id.csv".

So now you have the right order and the positions. What is stopping you from extracting the sequences?

By using the above mentioned command I formatted the coordinates in right order and made it compatible for samtools. But still, one more thing I have to do is reverse complement of those right ordered sequences. I have jumbled coordinates, so the command changes the order wherever it needed.

I need to extract the sequence based on the coordinates specified above (even though it is mixed up with inverse ordered and right ordered). But, Now I need to reverse complement the right ordered (right ordered by using above command) sequences. So that I could go for my further analysis. Moreover, I need to automate this process. The reverse complement implementation to be taken place, wherever the inverse coordinates are there.

What do you know biopython? It was designed to make reverse-complement sequences and translate if needed

The exonerate tool comes with some very handy utilities for sequence manipulation. Here, we are interested in fastasubseq and fasterevcomp. First, install exonerate as described here. If you're using Debian or Ubuntu or one of their derivatives, it might be enough to run
sudo apt install exonerate

Then, follow these steps:

Create one fasta file for each sequence in genome.fasta
fastaexplode genome.fasta

A little preprocessing
$ awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt
Contig3 7 8 -
Contig2 5 5
Contig1 3 9 -

So, that will print the target contig's name, the start coordinate of the target sequence, the length of the target sequence and a - if it is on the reverse strand. The fastasubseq program I will use below starts counting from 0 and not from 1, so we also need to subtract 1 from the coordinates to get them right.
Now, feed that into fastafetch and fastarevcomp
tmpFile=$(mktemp);
    > subseqs.fa
    awk -F[:-] '{
                if($2>$3){
                    start=$3-1;
                    len=$2-start" -"
                }
                else{
                    start=$2-1;
                    len=$3-start
                }
                print $1,start,len}' id.txt |
    while read cont start len rev; do
        fastasubseq "$cont".fa $start $len > $tmpFile;
        if [[ -n $rev ]]; then
            fastarevcomp "$tmpFile" >> subseqs.fa;
        else
            cat "$tmpFile" >> subseqs.fa;
    fi
    done

This will produce the following file when run on your data:
$ cat subseqs.fa
>Contig3:subseq(6,9) [revcomp]
GATCAGCGC
>Contig2:subseq(4,6)
AGCGCT
>Contig1:subseq(2,10) [revcomp]
CTTGATGGCC

Finally, remove the per-contig fasta files
rm Contig*fa

Thank you so much terdon for your time and help. Your script works awesome. However, I did a mistake in my question. Instead of mentioning Contig3:15-7, Contig1:12-3, I have mentioned Contig3:15-7 and Contig2:5-10.
I need to reverse complement the sequences, only if it in reverse order/inverse order. Which ever sequences are not in reverse order (Properly order sequences), only to be extracted as its not to be reverse complemented.

Please accept my sincere apologies for the same. Could you please do the needful.

@K.Dineshkumar I'm sorry, but I don't understand the problem. The script will only reverse complement if the start position is smaller than the end position. Isn't that what you want? In the example above, Contig3:15-7 and Contig1:12-3 have been reverse complimented, but Contig2:5-10 has not.

@K.Dineshkumar please edit your question and add this. Show us the expected output. It's very hard to read in the comments.

@K.Dineshkumar
Never mind, I found the problem. The fastasubseq tool uses 0-based coordinates (the first position of a sequence is 0), but you have 1-based (the first position of a sequence is 1). Try the updated answer.

@terdon Please see my edited question, wherein I have included my expected outcome.

@terdon, The revised/updated script printing one base less in the end of the sequence. For example,  Contig2:5-10 says, print the sequence from 5 to 10 But the script print only 5-9.

@K.Dineshkumar argh! Yes, of course, sorry. Fixed now, try the (last, I hope) update. And thanks for adding the output to your question, now I can check :)

@terdon, Now it works perfect. I really appreciate your help. Thank you.

@K.Dineshkumar you're welcome. If this answers solves your issue, please take a moment to accept it by clicking on the checkmark on the left. That will mark the question as answered and is the way that thanks are conveyed on the Stack Exchange sites. You don't need to do it now, you might want to wait in case anyone can give a better solution, but please remember to accept an answer since that marks the question as answered.

@K.Dineshkumar thanks for accepting! Feel free to accept one of the other answers though, I won't take it the wrong way :) They give more efficient solutions.

Since you tagged python, here is a solution that solves this in a single script using Biopython:
import sys
from Bio import SeqIO

d = {}
with open('id.txt') as f:
    for line in f:
        contig, coords = line.split(':')
        start, end = map(int, coords.split('-'))
        d[contig] = start, end

for record in SeqIO.parse('genome.fasta', 'fasta'):
    start, end = d[record.id]
    if start > end:
        record.seq = record.seq[end - 1: start].reverse_complement()
    else:
        record.seq = record.seq[start - 1: end]
    SeqIO.write(record, sys.stdout, 'fasta')

Output on your example:
>Contig1
CTTGATGGCC
>Contig2
AGCGCT
>Contig3
GATCAGCGC

Thank you @Chris_rands

I try to use it for a bacterial genome. But it shows error as follows,
Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.id] KeyError: '1'

@K.Dineshkumar try using d[record.description] instead, i can't know what will work without seeing your data though

My genome of interest can be obtained from the following link https://www.ncbi.nlm.nih.gov/Traces/wgs/CAHO01?display=contigs

My coordinate is as follows,                                                                                CAHO01000012.1:92353-95295 CAHO01000017.1:27362-27949 CAHO01000014.1:11449-12012 CAHO01000018.1:51156-52928 CAHO01000030.1:58079-56889 CAHO01000001.1:60348-62039 CAHO01000026.1:45206-44508

When I tried with d[record.description] instead of d[record.id], I end up the error as follows,                                                                                Traceback (most recent call last):   File "/home/ga/Documents/commands/scripts/amphora/reverse_comple_seq_extrac.py", line 12, in <module>     start, end = d[record.description] KeyError: 'CAHO01000001.1 Xanthomonas citri pv. mangiferaeindicae LMG 941 WGS project CAHO00000000 data, contig XMINcontig00001, whole genome shotgun sequence'

If you're okay with using a Perl script, you could convert your headers into stranded BED6 and use the following script bed2faidxsta.pl with samtools-indexed FASTA to get back stranded sequence:
https://gist.github.com/alexpreynolds/fa9b0f90e181e3b4f640
This script looks at the sixth column of your BED file of chromosome, start, and stop positions, to decide whether to do a reverse complement of the queried sequence.
You could use awk to quickly build such a BED6 file:
$ awk -vFS='[:-]' '($0~/^>/){ if ($2 < $3) { strand = "+"; } else if ($2 > $3) { strand = "-"; } else { printf("Coordinates must differ\n"); exit 1; } printf("%s\t%d\t%d\t.\t0\t%c\n", $1, $2, $3, strand; }' in.fa > out.bed6

To set up indexed FASTA, here's a generic shell-based approach, which uses [Kent utilities][1] fetchChromSizes to get a list of chromosomes for the build mm9:
$ mkdir mm9
$ cd mm9
$ while read chr; do echo "$chr"; wget -qO- "http://hgdownload.cse.ucsc.edu/goldenpath/mm9/chromosomes/$chr.fa.gz" | gunzip -c - > $chr.fa; samtools faidx $chr.fa; done < <(fetchChromSizes mm9 | cut -f1)

This can be adjusted for mm10, hg19, hg38, and so on.
Once you have indexed FASTA and the above Perl script, you could do something as easy as the following:
$ ./bed2faidxsta.pl --fastaDir="/path/to/indexed/mm9" < in.bed > out.fa

If you're working with a reference genome, the nice thing about this approach is that awk, Perl, and indexed lookups are fast. You can reuse indexed FASTA for other projects that also require sequence lookups by genomic interval.

Thank you @Alex Reynolds


text:
Discordance of Gene Ontology data between Mouse Genome Informatics and Ensembl sources
In my first approach to the use of Gene Ontology, I have found a discordance if we consult either Mouse Genome Informatics (MGI) or Ensembl for these data.
Concretely, for the mouse genes Zic1 (NCBI 2271) and Zic2 (NCBI 2272), MGI assigns them  up to 23 and 24 different GO terms, respectively.
However these numbers are different according to Ensembl data retrieved through R/Bioconductor. In particular, for Zic2 (Id 22772) the difference is very significant (24 terms in MGI vs. only 3 terms in Ensembl):
> my_ensembl = useMart ("ensembl", dataset = "mmusculus_gene_ensembl")
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22771', mart = my_ensembl)
        go_id
1  GO:0003676
2  GO:0046872
3  GO:0003677
4  GO:0005634
../..
22 GO:0008589
23 GO:0007389
24 GO:0021510
25 GO:0042472
> getBM (attributes = "go_id", filters ="entrezgene" , values = '22772', mart = my_ensembl)
       go_id
1 GO:0003676
2 GO:0005634
3 GO:0016604
>

So, my questions are: Which is the reason for this discordance? Has it something to do with the different types of evidence to assign GO terms? Also, in order to proceed to analysis like GO enrichment, which of these two databases would be more appropriate?

I'm looking into this. We should have all the annotations from MGI and I'm not sure why we don't.

Short answer, we get our annotations from GOA, which has the three listed. We're in touch with GOA to see why they differ from MGI.


text:
Getting a stretch of genomic ranges from a dataframe/granges object based on metadata column
I have a "test" data.frame object in R, which is basically a small subset of a 66000 row dataframe, which looks as follows:
chr   start    end      PC1     comp
chr1  3360001  3400000  -0.009  B
chr1  3400001  3440000  -0.004  B
chr1  3440001  3480000  -0.001  B
chr1  3480001  3520000  0.003  A
chr1  3520001  3560000  0.002  A
chr1  3560001  3600000  -0.001  B
chr1  3600001  3640000  0.002  A
chr1  3640001  3680000  0.004  A
chr1  3720001  3760000  0.003  A
chr1  3760001  3800000  0.003  A
chr1  3800001  3840000  0.007  A
chr1  3840001  3880000  0.006  A
chr1  3880001  3920000  -0.004  B
chr1  3920001  3960000  -0.017  B

As you can see, the interval between each start and end position is 40kb. The column comp can have only two values, A or B. Going from row 1 to the last row of my dataframe and based on the values of the column comp, I want to get the whole stretch/block of the genomic region, once the compartment(character in column comp) is changed, and perform a mean of the PC1 values in that stretch.
So in this example, I basically want to shrink my test dataframe as follows:
chr1  3360001  3480000  mean(c(-0.009,-0.004,-0.001))  B
chr1  3480001  3560000  mean(c(0.003,0.002))  A
chr1  3560001  3600000  mean(c(-0.001))  B  <---- #notice this is a 'singleton'
chr1  3600001  3880000  mean(c(0.002,0.004,0.003,0.003,0.007,0.006))  A
chr1  3880001  3960000  mean(c(-0.004,-0.017))  B

I tried with a simple code, but it only runs until the first compartment value change(so row 4 in this test), and I am unable to figure out how to make it run for the whole dataset of 66000 rows. Also the use of multiple for loops, which is what I tried, may not be the prettiest way to do this.
I guess there exists a Bioconductor package which performs this task very easily. If anyone can point me to that direction, it'd be great.

The GenomicRanges package in bioconductor will do this.  It lets you do 'set' operations on genomic ranges, such as concatenating them, finding overlaps, etc.  Hopefully something in the following code will be inspiring, even if I haven't solved the exact problem for you here:
library(GenomicRanges)

# make a GenomicRanges set corresponding to the original table
test.gr = makeGRangesFromDataFrame(df = test, start.field = "start", end.field = "end", seqnames.field = "chr")

# make separate GenomicRanges sets for each of the two partitions of the genome
test_A.gr = makeGRangesFromDataFrame(df = test[test$comp=="A",], start.field = "start", end.field = "end", seqnames.field = "chr")
test_B.gr = makeGRangesFromDataFrame(df = test[test$comp=="B",], start.field = "start", end.field = "end", seqnames.field = "chr")

# apply the reduce() function to concatenate segments which are adjoining each other into single segments
reduced_test_A.gr = reduce(test_A.gr)
reduced_test_B.gr = reduce(test_B.gr)

# group PC1 values by segment classification
test[queryHits(findOverlaps(test.gr, reduced_test_A.gr)), "PC1"]
test[queryHits(findOverlaps(test.gr, reduced_test_B.gr)), "PC1"]

Take a look at the rle function:

Compute the lengths and values of runs of equal values in a vector – or the reverse operation.

In this way you can mark the moment in which A changes to B and the length of the streches of A and B.
With that I think you can reduce the number of for loops you have to use (I can not think a way of doing this without at least one loop, but for sure it may exist)

What is proposed by @Jonathan works, just adding another possible solution here.
library(GenomicRanges)

test.gr = makeGRangesFromDataFrame(test,keep.extra.columns=TRUE)
# as long as the previous comp does not agree with current, we make it a new segment
test.gr$seg = c(0,cumsum(diff(as.numeric(test.gr$comp)) != 0))

# simple function to collapse the range
fn = function(u){
  out = reduce(u)
  out$comp = unique(u$comp)
  out$mean = mean(u$PC1)
  out
}

reduced_range = unlist(GRangesList(lapply(split(test.gr,test.gr$seg),fn))))

```


text:
Get residue position from Pymol Selection
I can make a selection in Pymol using the resi parameter. Is there a way to get this resi data back out of a selection object? The only thing I found so far is how to get the position, as in the x/y/z coordinates from a selection.
For context, I have a table with a set of data points for each residue in the amino acid sequence of the displayed structure and I want to write a function where I can print the data based on the current selection. So my thought would be something like the following, where the content of the 'get_residue_positions_for_selection' function would be what I'm asking about.
def show_selection_position_data():
    residue_data = {
        "header": ['variation', 'rmsf']
        , "1": ['high', 0.1]
        , ...
    }
    selections = cmd.get_names('selection')

    if sele not in selections:
        print("No selection")
        return

    residues = get_residue_positions_for_selection("sele")

    for residue in residue_selections:
        print(residue_data[residue])

    return

This certainly a very interesting question and a cool answer. The advice I was given (to pretty much the same problem) was to use Pymol scripting directly rather than via Python.

Ok, after some more searching I found the answer.
You have to use the cmd.iterate function to iterate over the default selection ("sele", first parameter). With the second parameter you can push the variables you need into an array for later access. In this case I wanted the residue number so I used the variable 'resv'. A full list of the variables exposed by the iterate function can be found on the pymolwiki
Note that you can use a function call as the second parameter and if you need a more complex selection than 'sele' you can use selector.process in the iterate function.
Here is my final working function. Note that I used a set for storing the residue ids, since the iterate function technically iterates over the atoms and so the same residue will likely appear multiple times in your selection.
def show_selection_position_data():
    residue_data = {
            "header": ['variation', 'rmsf']
            , 1: ['high', 0.1]
            , 2: ['high', 0.2]
            , 3: ['low', 0.1]
            , 4: ['high', 0.05]
            , 5: ['middle', 0.2]
            , 6: ['high', 0.1]
            , 7: ['high', 0.1]
        }
    selections = cmd.get_names('selections')

    if "sele" not in selections:
        print("No selection")
        return

    userSelection = "sele"

    stored.residues = set()
    cmd.iterate(selector.process(userSelection), 'stored.residues.add(resv)')

    for residue in stored.residues:
        if residue in residue_data:
            print(residue_data[residue])
        else:
            print("{}: no data".format(residue))

    return

Nice. And thank you for taking the time to post an answer!


text:
Why does DESeq2 convert numeric columns to factor during differential expression analysis?
I'm attempting to perform differential expression analysis using DESeq2 on a dataset of ~90 samples and ~20000 transcripts. I have a number of variables I'd like to test against, some of which are qualitative (factors) and some of which are quantitative (numeric). A sample of my data are below:
> head(data.phenotype.test)
# A tibble: 6 x 7
  Novogene_ID BORIS_ID     treatment queenness ovaries   elo sampletime
  <fct>       <fct>        <fct>         <dbl>   <dbl> <dbl> <fct>
1 W03BLGR     003W-blu-grn Qr12        0.0202    124.  1000  12
2 W03GR       003W-grn     Qr12        0.00210    66.8  923. 12
3 W03SL       003W-sil     Qr12        0.969     316.  1077. 12
4 W04BLOR     004W-blu-ora Qr12        0.00168    66.4  876. 12
5 W04GRSL     004W-grn-sil Qr12        0.175     210.   931. 12
6 W04YLWH     004W-yel-whi Qr12        0.984     302.  1192. 12

Now, let's say I want to test which transcripts are differentially expressed together with a numeric variable such as 'ovaries':
dds.transcript.ovaries = DESeqDataSetFromMatrix(countData = data.transcript.count.ovaries,
                                        colData = as.matrix(data.phenotype.ovaries),
                                        design = ~ ovaries)

This should work, but DESeq2 does something strange- it converts 'ovaries' to a factor, with the message "some variables in design formula are characters, converting to factors". But there are no character values in the supplied data frame, and certainly not in the column 'ovaries'! This conversion to factors is a big problem for me because, having converted a numeric column to factors, DESeq2 unsurprisingly reports that differential expression analysis cannot be performed because there are an equal number of factor levels to samples.
I have found a temporary workaround: if I reduce the data frame to just the 'ovaries' column, DESeq2 no longer converts the numeric data to factor levels and I'm able to perform differential expression analysis as normal. The problem is that I'd eventually like to perform multivariate analyses, e.g. with a design like:
~ ovaries + elo + treatment

However, doing so isn't going to be possible if DESeq2 keeps converting the numeric columns to factors, which it persistently does unless the numeric columns are supplied alone (which of course they can't be for a multivariate design).
What's going on here? Why does DESeq2 convert the numeric columns to factors, seemingly having misinterpreted those columns as characters?
NB: I have replicated this issue with different data on a different machine and the result is the same, so I don't believe this an issue with my R setup.
Edit: Here's the structure of the df:
> str(data.phenotype.test)
Classes ‘tbl_df’, ‘tbl’ and 'data.frame':   64 obs. of  7 variables:
 $ Novogene_ID: Factor w/ 64 levels "W03BLGR","W03GR",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ BORIS_ID   : Factor w/ 64 levels "003W-blu-grn",..: 1 2 3 4 5 6 7 8 9 10 ...
 $ treatment  : Factor w/ 4 levels "Qc","Qr12","Qr3",..: 2 2 2 2 2 2 2 2 2 2 ...
 $ queenness  : num  0.02017 0.0021 0.96853 0.00168 0.17529 ...
 $ ovaries    : num  124.1 66.8 316.1 66.4 209.9 ...
 $ elo        : num  1000 923 1077 876 931 ...
 $ sampletime : Factor w/ 3 levels "0","3","12": 3 3 3 3 3 3 3 3 3 3 ...
 - attr(*, "na.action")= 'omit' Named int  1 2 3 4 5 6 7 8 9 10 ...
  ..- attr(*, "names")= chr  "1" "2" "3" "4" ...

Edited:
It's the "as.matrix" that is causing the problem.  When I use that on my sample annotation data, and make dds with that, it turns my numbers into factors

Hi, thanks for your reply. The actual values in the ovaries column do continue past the decimal point, they're just being cut off in the tibble display. I'll add the results of str(data.phenotype test) to the OP.

You were right! Removing the as.matrix part solved the problem. Thank you!


text:
Generating DNA sequences with constraints
I would like some advice on potential strategies to address the following problem.
I want to write a program that will generate DNA sequences that are optimized on two constraints based on an input genome:

The GC content should be as close to the genome as possible.
When taking all k-mers (of a predefined length) in the sequence, the sum of their frequencies in the genome should be as small as possible (i.e. it should be made of rare k-mers).

I currently score sequences using  :
$S = w_{GC} * (1 - |GC_{genome} - GC_{seq}|) + w_{Kmer} * (1 - \frac{mean(f(Kmer_{seq})}{max(f(Kmer_{genome}))}$
Where $w_{GC}$ and $w_{Kmer}$ are user-defined weights defining the relative importance of GC content and k-mer profile) and $f(Kmer_{seq})$ is the frequency of k-mers in the sequence.
Of course I could just use bruteforce and generate many random sequences of appropriate GC content, score them all and pick the very best. But 1) that's no fun and 2) it's slow.
The first thing I attempted was to extend the sequence one basepair at a time. Each iteration has 4 candidates (with a new A, C, T or G added). I pick the candidate using their score ($S$) as probability weight. This "blind" method is too simple and does not yield satisfactory results.
I thought a dynamic approach could work well (filling a table and backtracking),  but I am not sure this is applicable since the cost of extending the sequence would be dependent on all new k-mers generated and not be a fixed value (such as in sequence alignment).
I also thought something like reinforcement learning could be appropriate here, but I don't know much about it, so I am not sure.
I know all of this is probably overkill, but this is mostly for the challenge / learning :)
Which strategy would be most appropriate for this ? Did I miss something simple ? All advice / suggestions are welcome !

What is the purpose of simulating sequences from a genome, rather than simply using subsequences from the genome? If this is purely a programming exercise, the question would be better asked on a different StackExchange site. If there is a biological investigation behind it, you may be ignoring genomic context that is necessary for the proper function of the DNA sequence (e.g. see here).

@gringer I want to generate sequences that will minimize homology to the target genome at the scale of the kmer, so I don't want to just pick subsequences. Right, it might have been a good question for stack overflow as well !

You can use a $l$-order Markov chain. Here is the procedure:

Count $l$-mers in your genome. For small genomes, you can do that in Python. For large genomes, you may need jellyfish or KMC3.
Draw a $l$-mer randomly based on the distribution of $l$-mers. This generates the first $l$ bases.
Let $s$ be the last $(l-1)$ subsequence from the generated sequence. There will be at most four $l$-mers prefixed with $s$. It gives you a distribution of the next A/C/G/T base conditional on $s$. Randomly draw a base from this distribution. You get the next base.
Repeat step 3 until you generate a sequence long enough.

The generated sequence will have the same $l$-mer distribution to the original sequence. The GC content will be the same, too. In your case, if $k$ is small enough, you can let $l=k$. If $k$ is large, you have to choose a smaller $l$.
The procedure above doesn't consider the strand symmetry of DNA, but probably this doesn't matter much. It is also possible to adapt the algorithm to address this.
As to implementation, you can use python, but a high performance language like C will be tens of times faster at least.


text:
Deep learning RNA sequences
Currently I'm working on a project, which combines deep learning with RNA sequences. I'll try to predict pseudotorsion angles [1] from raw rna sequence. The ideas is to train a neural network with raw rna sequences and for each nucleotide their corresponding pseudotorsion angles, and than to predict the angles of the remaining sequences in the test set.
How my data is structured:
example:
seq1:   A  C  G    G   U   A   C
Eta:   169 87 110 87  45  187 78
Theta: 123 10 45  168 132 34  100

[1] These are angles describing the backbone conformation of a rna molecule.
I'm pretty new to the field of deep learning, and so far I build a simple feedforward neural network, but it's prediction accuracy is pretty low with only one percent.
Has anyone some tips for me how to improve this?
How do I preprocess this kind of data correctly for deep learning?
I appreciate any help.

Note that RNA has over 110 known modifications, most of which can't be detected using existing automated sequencing tools. These would probably influence the structure of the molecule, and reduce the accuracy of your prediction.

The question of folding of RNA sequences looks slightly similar to protein folding - perhaps searching in this domain might bring more suggestions.
An example of (current) state of the art of deep learning approach to protein folding is AlphaFold developed by DeepMind for CASP competition . They basically decompose the angles to an image representation and look for structures in the proximity of each atom. While it requires a significant processing power, RNA structures might be easier due to the size of the sequence.

I come from a protein background and this problem is analogous to protein torsion angle prediction, which in turn is a variant of protein secondary structure prediction.
Conventional ways to go here would be to use a sliding window of x bases either side to make a prediction for each point, or use a CNN/LSTM-type model to run over the whole sequence and output angles.
There are other technical concerns such as how to split the training and test sets, how to output the angles from the model (usually transformed by sin or cos), and whether to predict a single value or the probability of the value being in certain ranges.
The reference to protein tertiary prediction and CASP in another answer isn't directly relevant here unless you intend to build 3D models from your torsion angles, which you might very well do.
Some references on torsion angle prediction for proteins:

A recent review.
A recent state of the art method.
A classic paper describing some of the approaches.

Thanks for your post and highlighting your work. Could you clarify that CNN is convolution neural network and when do you use it as opposed to LHTM?

Correct, CNN is convolutional neural network and LSTM is long short-term memory, a variant of recurrent neural networks (RNNs). LSTMs would typically be used for sequential data, where you aim to predict a property for each point in the sequence or predict the next point in the sequence.
CNNs apply the same filters to every point in the sequence and rely more on picking up local features than treating the input as a sequence. They are usually applied to images (2D) but can be 1D, 2D, or 3D (etc.) as required.

LSTM, even I know this acroynm. Many thanks that finally makes  sense in context to how CNNs operate in this field and more generally. We're using your (DJ's) stuff, it always nice to know how it works :) . For the OP getting deep learning to work isn't always trivial, this team is expert.

I think the solution is more in the DL.
I would look at simpler ML, using a random forest as a starting point and then maybe SVC. If you get some reasonable predictive power on this then move towards neural network/Tensorflow type calculation.


text:
Error Code in WGCNA: Dealing with large data sets: block-wise network construction andmodule detection
I am comparing single block module detection with 2 block module detection in WGCNA.
Does anyone know what this relates to? I have rechecked and rechecked my codes for both the single and 2 module detection and cannot see an error. Thank you
I get the error code:
Error in plotOrderedColors(dendro$order, colors = colors, rowLabels = rowLabels,  :
  ERROR: length of colors vector not compatible with number of objects in 'order'.

Using the following code:
net= blockwiseModules(Th1_h5, power=6,
                      TOMType = "signed", minModuleSize=30,
                      reassignThreshold =0, mergeCutHeight = 0.25,
                      numericLabels = TRUE, pamRespectsDendro = FALSE,
                      saveTOMs = TRUE,
                      saveTOMFileBase = "PBCreponseTOMB",
                      verbose= 3
mergedColors = labels2colors(net$colors)

plotDendroAndColors(net$dendrograms[[1]], mergedColors[net$blockGenes[[1]]],
                    "Module colors",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang = 0.05)
moduleLabels = net$colors
moduleColors = labels2colors(net$colors)
MEs = net$MEs
geneTree = net$dendrograms[[1]]

bwet = blockwiseModules(Th1_h5, maxBlockSize = 8000,
                        power=6, TOMType = "signed", minModuleSize = 30,
                        reassignThreshold = 0, mergeCutHeight = 0.25,
                        numericLabels = TRUE,
                        saveTOMs = TRUE,
                        saveTOMFileBase = "Th1_NR_R_TOM-blockwise",
                        verbose=3)
bwLabels = matchLabels(bwet$colors,moduleLabels)

bwModuleColors = labels2colors(bwLabels)

plotDendroAndColors(geneTree,
                    cbind(moduleColors, bwModuleColors),
                    c("Single block", "2 blocks"),
                    main = "Single block gene dendrogram and module colours",
                    dendroLabels = FALSE, hang = 0.03,
                    addGuide = TRUE, guideHang =  0.05)

You will need to provide more detail. How big was the data from which you calculated the network (how many genes) and how did you determine geneTree?

I have 12711 genes after filtering and have added the additional code above. It has worked with the simple Automatic, one-step network construction and module detection and Step-by-step network construction and module detection but not when trying with 2 blocks. Thank you kindly

The problem is that blockwiseModules split your data into (probably 3) blocks, because the default maxBlockSize is 5000 and smaller than the number of genes in your data. If you have enough RAM (16GB is recommended for 12k genes, but 8GB just might be enough), rerun blockswiseModules again with maxBlockSize raised to above the number of your genes. Since you were able to run the function with maxBlockSize=8000, hoepfully you'll be able to do the single block calculation.

Thank you so much! I will give it a try now, thank you!

It works! Thank you so much, as always for your time and help!


text:
Issue with visualising cladogram/phylogenetic tree with multiple sequence alignment data in R?
I would like to visualize tree with multiple sequence alignment.
My din.newick file is shown below,
(org1:0.49999997,org2:0.47298786,org3:28.37858179);

My din.fasta alignment file shown below,
>org1
ATGCTA
>org2
TAGCAT
>org3
AGTAGT

I have tried the following code,
library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta")

The above mentioned code is generating tree with the alignment matrix, but the tree root tip labels are not visible. It is either overlapped with matrix or it is not at all visible.  Kindly help me to fix this issue.

You were very close, you just need to supply the offset argument to msaplot().
So using your code block above with one extra argument in msaplot()
library(tidyverse)
library(ggtree)
library(seqinr)
tree = read.tree("din.newick")
msaplot(p=ggtree(tree), fasta = "din.fasta")
msaplot(p=ggtree(tree) + geom_tiplab(align=TRUE), fasta = "din.fasta", offset = 5)

yields

Thank you for supplying a minimal and reproducible example!

Thank you @Kohl Kinning.

I suspect the solution is to realign the two plots within "ggarrange", available here
The problem is stacking ggtree with geom_tiplab. The ggarrange library is described on STHDA here
I have not used this library however, but it should be easy because you are working singly within R. I tend to use ggplot2 within Python via rpy2, which is much more difficult and restricted.


text:
Occupancy of TFs with the target genes
The occupancy of SMARCD3 in the target genes listed below. I want to see average, normalized ChIP-seq signal at the promoter proximal region (1000bp upstream and downstream of the TSS).
I have 4 different experimental conditions (overlayed in one plot, per target gene using different color codes) where we want to visualize the changes in occupancy of GATA1 under all the conditions at individual gene promoters. The following are our target genes

ABCB1
ABCC1
ABCC2
SNAI2

So how GATA1 occupancy is changed across all these target genes across 4 different treatment condition.

So far what i have done

Aligned the samples
Done peak calling
Annotated the peak

I did steps 2 and 3 using homer.
Now after all this how do I find the occupancy of a TF with its target genes? Do I get the information from my peaks annotated output or from the peak files?

I would ignore peak calling for this and instead compute enrichment of ChIP/input for the genome (e.g., with deepTools or presumably homer) and then plot it for the genes of interest individually (e.g., using IGV or pyGenomeTracks) or as a group (e.g., with computeMatrix). If the peaks are obvious and you trust your peak calling then sure you can just use that instead, but given that that's prone to error it's usually better to simply use normalized signal.
Regarding actually looking at occupancy, I don't think that's useful. Occupancy isn't binary. If you're instead wanting to associate peaks with genes then using the nearest gene is common (it's not great, but it's not terrible either).

actually i looked into your deeptool as many publications used to similar analysis but not sure how it was done.  " Occupancy isn't binary. " yes i now get it occupancy can;t be binary , " If you're instead wanting to associate peaks with genes then using the nearest gene is common" how do i perform this part in terms of tool can i  do it in deeptools?

For finding the nearest gene, just use bedtools nearest from bedtools.

https://www.sciencedirect.com/science/article/pii/S109727651930228X in this paper this figure (G) Metaplots of canonical PRC1 (PCGF2, CBX7, and PHC1) cChIP-seq at classical Polycomb chromatin domains (n = 2,096) in Pcgf4−/−;Pcgf2fl/fl ESCs (UNT and OHT), in this paper did they get into this? did they use bedtools nearest as well to find the gene as they write "recruits canonical PRC1 to Polycomb target genes in order to compact chromatin "

They didn't do any gene assignments in that subpanel, they just input known polycomb domains.

okay can you tell me how did they get the profile of particular genes in other words where do get this information for doing Metaplots? is it from the peak annotation

Those are from igv or pygenometracks or similar.

okay may be the data i have used the read counts are really low so i dont get that profile

I hope you are talking about this one G) Metaplots of canonical PRC1 (PCGF2, CBX7, and PHC1) cChIP-seq at classical Polycomb chromatin domains (n = 2,096) in Pcgf4−/−;Pcgf2fl/fl ESCs (UNT and OHT).

No, they didn’t need or use genes for subpanel G.


text:
Find intervals that genes fall within their range
I am trying to run some frequency based stats to identify selective sweep in my study system using Rpackage "PopGenome".
To proceed with it I have split my genome data (whole genome sequencing) into 5MB chunks
cregions <- character(31)
cregions[1] <- "1-5000000"
cregions[2] <- "5000001-10000000"
cregions[3] <- "10000001-15000000"
cregions[4] <- "15000001-20000000"
cregions[5] <- "20000001-25000000"
cregions[6] <- "25000001-30000000"
cregions[7] <- "30000001-35000000"
....

I have a bunch of genes and I want to know each of my genes falls within the range of which intervals in the above list.
these are my genes:
chrom start stop gene_id
Chr01 107310787 107323633 gene2748
Chr01 113789595 113794851 gene2894
Chr01 159041389 159068535 gene4806
Chr02 2742629 2758060 gene4846

I'd like to get an out put like this
gene2748  window_2
gene2894  window_3
...

Any idea how I can do it in R?

The GenomicRanges package is a great place to start, something like this:
library(GenomicRanges)

# make a GenomicRanges set corresponding to the genes
genes.gr = makeGRangesFromDataFrame(df = genes, start.field = "start", end.field = "stop", seqnames.field = "chrom")

# turn your windows into a data frame and make a GenomicRanges set from that
windows = ...
windows.gr = makeGRangesFromDataFrame(df = windows ...

# find overlaps between the two sets
olaps = findOverlaps(genes.gr, windows.gr)

# output relevant fields
cat(paste(genes[queryHits(olaps),"gene_id"], windows[subjectHits(olaps),"window_id"]))


text:
What is the best way to start with Structural bioinformatics with zero prior experience?
My undergraduate thesis is on the topic "Ligand design from protein-protein interfaces" and my PI has asked me to learn the basics of Structural bioinformatics over the summer so that I'll be able to start with my project as I begin my next semester(last week of July). Can someone please suggest me top resources to learn the basics of structural bioinformatics. Though I have done projects in Genomics, and Systems biology, this will be my first time with structural bioinformatics, so, I'll really appreciate some good leads. Thanks!

Hello Akhil, this is not the greatest place for open questions. You could try one of the forums or just google it.

As others have said, there is no short cut and you will need to spend some time with a good book.
I would say there are 3 broad areas of understanding you have to have as a structural bioinformatician:

Solid understanding of the basic science. This includes molecular biology but also chemistry, which is often overlooked.
Knowledge of the databases and data formats available in the field. You need to know where to go to get protein structures, protein sequences etc.
Programming and algorithmic skills. You don't have to be a computer whiz but being effective in a high-level language (e.g. Python, R, Julia) and being able to parse data and carry out basic data analyses are essential.

If you want to take it further and have a bit more time then I would consider implementing something yourself such as a protein secondary structure predictor. Getting this working is challenging and would show you have competency in the above three areas.

Thank you. If you could also suggest me some resources(your top picks) for the first two areas, it'll be really helpful.

Of the top of my head "Introduction to Protein Structure" by Branden and Tooze is good. In terms of databases have a look at PDB, UniProt, Pfam, CATH and InterPro.

Thanks @jgreener

A good starting point could be this book by Gu and Bourne:
Structural Bioinformatics, 2nd Edition

Thanks, I'll check it out. But I was hoping for a good blog or some tutorial(web/video) kind of resource rather than a book

@AkhilVerma — Books are far more informative. They contain more information, undergo quality control if published by reputable publishers and are generally better written, more comprehensive and more reliable. You cannot learn anything worthwhile without putting in some effort. Certainly not structural biology. How’s your chemistry, incidentally? You’re going to need that too


text:
Function to find ATP binding sites
I'm currently doing a python project where I need to write code that extracts the ATP binding site from a FASTA file using a regular expression.
I've tried the following code but it does not work. It opens a new file but the file is empty. I don't know what's wrong.
from Bio import SeqIO
import re
import os
path="C:/Users/nessr/Desktop/info projet"
os.chdir(path)
records=list(SeqIO.parse("plasmo.fasta","fasta"))
ATPM="[LIV]G[^P]G[^P][FYWMGSTNH][SGA][^PW][LIVCAT][^PD].[GSTACLIVMFY].{5,18}[LIVMFYWCSTAR][AIVP][LIVMFAGCKR]K"
patern=re.compile(r"(ATPM)")

def atp_binding(records,patern):

    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()

            iend=iend+50

            print(istart,m.group())

file = open("output.txt","w")

for seq_record in SeqIO.parse(records, "fasta"):
    file.write(atp_binding(seq_record.records,patern)+"\n")
    file.write('-------------------\n')
file.close()

Welcome to Bioinformatics.SE! We're happy to help with homework and research problems, but you're much more likely to get helpful feedback if you describe precisely what your problem is. Asking "can you please tell me how to do this?" will rarely get you a useful answer. Dumping poorly formatted code is another way to turn potential helpers away. Please take a moment to clean up your Python code (use the Markdown Editing Help guide if needed), and then please be more specific about what is and is not working.

I put backticks around your code and cleaned up the question a little bit. I hope this helps. Please continue to edit and clarify your question.

The problem in your code is that the function atp_binding is not returning any value, but printing it to the stdout.
You should add a return statement with the result you want to save in your output file.
Also, you should change how the resulting value is being saved. As the return value of your function is None, the file handler will write nothing.
I would use something similar to this:
def atp_binding(records,patern):
    s = ""
    for item in records:
        for m in patern.findall(str(item)):
            istart=m.start()
            iend=m.end()
            iend=iend+50
            s+="{0}\t{1}\n".format(istart,m.group())
    return s

file = open("output.txt","w")

for seq_record in SeqIO.parse(records, "fasta"):
   file.write(atp_binding(seq_record.records,patern)+"\n")
   file.write('-------------------\n')
file.close()

It is a good practice to handle file with with command. like : with open("output.txt", "w") as file : .....


text:
Import gene list to Seurat to define cell types
Is there a way to import gene list into Seurat to define cell type? The default cell types in Seurat is not enough for our research. For example, we want to mark a subtype of B cells in Seurat, but seems like Seurat only have general B cell, I'm planning to download gene list from Genomic Cytometry, and then import the gene list into Seurat and then Seurat define cell types in the cluster based on the gene list I import. Does Seurat support this function?

Welcome to the site. Did you find anything on the vignette and on the several questions (and answers) posted online?

I checked, but unfortunately I couldn't find any information.

I am in roughly the same situation: I have a list of marker genes for each cell type, and I now I want to assign a cell type to each cell using this. Did you ended up finding solution and would you mind sharing?

You can use scID with your gene list to find the cells that match your gene list. You can then find which cluster in Seurat have these cells.

Seurat does not define cell types by name. It clusters and assigns each cell to a cluster, from 0 to X.
If your data has the cell type (e.g. B,T, Mast cells) it means that someone annotate the clusters so that they have a biological meaning.
You can assign different names to the clusters by using the AddMetaData function.
sampleID  cellType
cell1     B
cell2     B
cell3     T
cell4     T

newMeta <- oldMEta
newMeta$origType <- newMeta$cellType
newMeta$cellType <- c("B1","B","Treg","Tkill")

object <- AddMetaData(object, newMeta)

The new metadata should look like
sampleID  origType  cellType
cell1     B         B1
cell2     B         B
cell3     T         Treg
cell4     T         Tkill

You can then swap the labels using SetAllIdent(object = test, id = "cellType"). Your tSNE or uMAp, will now show the new cell types.
Have a look at Manually define clusters in Seurat and determine marker genes

Thank you so much for your answer. But actually I want to define a new subtype of B cells in the B cell cluster. I know the genes that express the new subtype of B cells, but the current Seurat cluster only showing the general B cell. I'm thinking about doing a subtype of B cells, and then import the the gene list to the Seurat to see if Seurat can show the subtype of B cells based on the gene list I provided. I'm not sure if Seurat can do this.

Seurat does not use gene lists to cluster the data. What I usually do to look for sub-cell types, is to subset the dataset (only the cells from one cluster) and re-run the clustering/PCA steps to see how those cells (and only those cells) are clustering together. The expression of your genes will tell you if the sub-clusters contains specific sub-types.  Have a look at the Seurat documentation under "Selecting particular cells and subsetting the Seurat object".

Hello, you mentioned that Seurat not use gene list to cluster data, but I thought that Seurat cluster cells based on the expression of certain genes, if a cells have high expression in certain genes, they will be define as a specific cell types, do you mean that this is not the case?

Seurat clusters based on the expression of certain genes, but those genes will be determined independently for each experiment. By default, it would be the most variable genes. It does not know if those genes have been previously described to be associated with a specific cell type.

Maybe you can try Seurat::AddModuleScore(), then FeaturePlot() and see if some of your B cells are different. After plotting this on GenePlot(), perhaps you can set a cutoff, then assign identities.
Alternatively, use your B cell gene list in RunPCA(object, pc.genes = yourgenelist) instead of the usual variable genes. You can also subset B cells first with SubsetData(), and run the above.


text:
Any way to quantify the variation of genes that expressed in Affymetrix expression data?
am experimenting preprocessed Affymetrix microarrays expression data matrix (Affymetrix probe-sets in rows (32830 probesets), and RNA samples in columns (735 samples)) for my downstream analysis. Here is how the data looks like:
> dim(eset_HTA20)
[1] 32830   735
## load assayData
HTA20_rma <- load("HTA20_RMA.RData")

> eset_HTA20[1:10,1:3]
             Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at                6.062215        6.125023        5.875502
10_at               3.796484        3.805305        3.450245
100_at              5.849338        6.191562        6.550525
1000_at             3.567779        3.452524        3.316134
10000_at            6.166815        5.678373        6.185059
100009613_at        4.443027        4.773199        4.393488
100009676_at        5.836522        6.143398        5.898364
10001_at            6.330018        5.601745        6.137984
10002_at            4.922339        4.711765        4.628124
10003_at            2.689344        2.771010        2.556756

objective:
I am interested in to reduce the dimension of the dataset by tossing off low expressed in genes in the experiment. To do so, I need to quantify the variation of the gene that expressed. I checked the possible approach and using the coefficient of variation (cv) could be one of the approaches I could try. It is not intuitive for me how to quantify the variation of genes that expressed.
Can anyone point me out how to lay out this procedure on the above expression data? Is there any Bioconductor package to do this kind of task? any idea to make this happen? thanks

You can use the following code to calculate the coefficient of variation:
# expr is your expression matrix.
SD <- apply(expr, 1, sd)
CV <- sqrt(exp(SD^2) - 1)

It might be implemented in some package but it is so brief that you can write it again yourself. Then you can filter out those that are below certain percentage of the distribution of CV (like 0.1)
expr[CV > quantile(CV, probs = 0.1), ]

Another way to reduce dimensionality of your dataset is via principal component analysis or PCA. It is implemented as pcaMethods R package, for instance.


text:
what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?
What is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets nowadays?
About 5 years or so ago, it used to be the latest release of the CSHL HapMap 12 populations that were part of the 1000 genomes project. These were:
CHB CHD MEX GIH TSI LWK ASW MKK HCB JPT CEU YRI
For example, the EnsEMBL project currently has these 12 populations available as a vcf here:
http://ftp.ensembl.org/pub/current_variation/vcf/homo_sapiens/

CSHL-HAPMAP-HAPMAP-CHB.vcf.gz
CSHL-HAPMAP-HAPMAP-CHD.vcf.gz
CSHL-HAPMAP-HAPMAP-MEX.vcf.gz
CSHL-HAPMAP-HAPMAP-GIH.vcf.gz
CSHL-HAPMAP-HAPMAP-TSI.vcf.gz
CSHL-HAPMAP-HAPMAP-LWK.vcf.gz
CSHL-HAPMAP-HAPMAP-ASW.vcf.gz
CSHL-HAPMAP-HAPMAP-MKK.vcf.gz
CSHL-HAPMAP-HapMap-HCB.vcf.gz
CSHL-HAPMAP-HapMap-JPT.vcf.gz
CSHL-HAPMAP-HapMap-CEU.vcf.gz
CSHL-HAPMAP-HapMap-YRI.vcf.gz

Each SNP has an AF entry from which a multi-populations vcf with rsids, alleles and frequencies can be built, where the AFs are such as AF_CHB, AF_CHD, etc.
The 1000 genomes project populations documentation describes more than these 12 populations, but I haven't seen equivalent population vcfs with AFs built from the individuals within the population for the remainder of these, apart from the original HapMap 12 marked with a * below:
###Populations and codes

*         CHB   Han Chinese             Han Chinese in Beijing, China
*         JPT   Japanese                Japanese in Tokyo, Japan
         CHS    Southern Han Chinese    Han Chinese South
         CDX    Dai Chinese             Chinese Dai in Xishuangbanna, China
         KHV    Kinh Vietnamese         Kinh in Ho Chi Minh City, Vietnam
*         CHD   Denver Chinese          Chinese in Denver, Colorado (pilot 3 only)

         CEU    CEPH                    Utah residents (CEPH) with Northern and Western European ancestry
*         TSI   Tuscan                  Toscani in Italia
         GBR    British                 British in England and Scotland
         FIN    Finnish                 Finnish in Finland
         IBS    Spanish                 Iberian populations in Spain

*         YRI   Yoruba                  Yoruba in Ibadan, Nigeria
*        LWK    Luhya                   Luhya in Webuye, Kenya
         GWD    Gambian                 Gambian in Western Division, The Gambia
         MSL    Mende                   Mende in Sierra Leone
         ESN    Esan                    Esan in Nigeria

*         ASW   African-American SW     African Ancestry in Southwest US
         ACB    African-Caribbean       African Caribbean in Barbados
         MXL    Mexican-American        Mexican Ancestry in Los Angeles, California
         PUR    Puerto Rican            Puerto Rican in Puerto Rico
         CLM    Colombian               Colombian in Medellin, Colombia
         PEL    Peruvian                Peruvian in Lima, Peru

*         GIH   Gujarati                Gujarati Indian in Houston, TX
         PJL    Punjabi                 Punjabi in Lahore, Pakistan
         BEB    Bengali                 Bengali in Bangladesh
         STU    Sri Lankan              Sri Lankan Tamil in the UK
         ITU    Indian                  Indian Telugu in the UK

I presume there is more public data nowadays that can be accessed to build a more complete vcf files of allele frequencies per population, with as many populations as possible to use as a reference dataset with new data.
Thanks in advance.

The best set I'm aware of can be downloaded from https://gnomad.broadinstitute.org .

Thanks, I have seen there is a Downloads section which hopefully will have the files that I am after.

I could only find continent-level AFs for gnomAD. I presume nobody built the population-level frequencies file for these individuals, but maybe someone did and this is available on a website?

Can you explain what you mean by "population-level?" I see ones for racial and other ethnic groups as well.


text:
How to use extractfeat?
I am working on a script on bash linux to get the CDS (coding sequence) of a gene using extractfeat by EMBOSS. It gives me the error: Warning: No sequences written to output file
I am unsure if my syntax is right, but this is what I have:
extractfeat inputFile outputFile -type CDS

Hi, KA97, what is the error? If you can include it, this will help establish what the problem is.

. The error is :   Warning: No sequences written to output file
It produces a blank output file with no sequence in.

Perhaps,removing introns (if it is a higher eukaryote)
extractfeat tembl:x65921 -type CDS -join

It is difficult to know without knowing the gene.
BTW NCBI's Entrez "efetch" will do essentially the same thing, but its more complicated

I tried that syntax with the -join at the end which didn't work either
If it matters, I am currently testing it on a Histone H1 gene of C. elegans
I'd prefer If I could get the CDS by processing the sequence file in my bash script. Since I'm using efetch to download the sequence anyway, I don't mind trying that out if it's not too complicated.

Okay, the syntax you are using is correct, but you don't need to specify the output file

I managed to fix it. I was accidentally using FASTA format files instead of Genbank flat file format which seemed to be the problem.

That's good: fasta format has no features for extractfeat to extract, so it would have no chance. The error trapping in the code is clearly absent and it might be worth reporting the bug to the Emboss team.


text:
What are the mathematical assumptions and algorithmic procedure behind MACS2?
I wonder if there's some mathematical details that outline what the algorithms does.

Hi and welcome, it would be helpful if you could put your question in context. What is the objective, how does the algorithm help and why it's published information is insufficient for you.

Have you read the papers?

See here and here.

@TheSimpliFire Thanks.

authors: ə̷̶̸͇̘̜́̍͗̂̄︣͟


text:
Biological replicates on Chip-seq Transcription factor data
I have 2 biological replicates for chip-seq transcription factor data and I have applied macs2 peak calling for each replicate separately.
How can i make use of the biological replicates to extract the most out the data?
Thanks!

Hello Dimtiris, that's quite a broad question, is there a way to narrow it? What authors of macs2 say? How it's treated in other publications? How do you think you should treat your two replicates?

You could first look at the degree of correlation between the two replicates - what proportion of peaks are shared between the two, versus peaks found in one sample only?  This will give an idea of how repeatable the analysis is, and how many peaks are variable between samples or due to artefacts of the method.
When it comes to interpreting the biology behind the binding patterns, you are probably going to want to only interpret peaks which you find in both replicates.

This question is somewhat generic, so a generic answer is that ENCODE has a Transcription Factor ChIP-seq Data Standards and Processing page that can give you a useful starting point.
For TF ChIP-seq data with replicates, the Irreproducible Discovery Rate (IDR) method helps leverage replicates to produce higher confidence peak calls, producing both "optimal" and "conservative" peaks, defined by some IDR threshold. Here's ENCODE's description of the method:

A statistical procedure that operates on the replicated peak set and compares consistency of ranks of these peaks in individual replicate/pseudoreplicate peak sets. Peaks with high rank consistency are retained. IDR can operate on peaks across a pair of true replicates resulting in a “conservative” output peak set, or across a pair of pseudoreplicates resulting in an “optimal“ output peak set. Peaks in the conservative peak set can be interpreted as high confidence peaks, representing reproducible events across true biological replicates and accounting for true biological and technical noise. Peaks in the optimal set can be interpreted as high-confidence peaks, representing reproducible events and accounting for read sampling noise. The optimal set is more sensitive, especially when one of the replicates has lower data quality than the other.

The Python code to run IDR on peak sets can be found on ENCODE's GitHub Chip-seq pipeline repo.


text:
Sibling vs Parent/Child relationship detection using whole-genome sequence data
I am analyzing a family-based whole-genome sequence (WGS) dataset.  My identity-by-descent(IBD) analysis shows they share about 50% of their DNA. I think my pedigree structure may have gotten mixed up though and I want to see if one individual is the mother or the sister of a male individual.  Is there anyway to do this either with Plink or KING?

With at least three samples, you can determine this by looking at the IBS0 output column (from either plink2 --make-king-table or KING): if the relationship in question has a very low IBS0 value compared to the other pairs, that's parent-child.  You can try to do the same thing with only two samples, but with no other IBS0 values to compare against it's easier to misjudge.


text:
Any way to filter out highly correlated genes with limma linear model?
I have Affymetrix gene level expression matrix (genes in the rows and sample ID on the columns), and I have annotation data of microarray experiment observation where sample ID in the rows and description identification on the columns.
example data
here is Affymetrix gene level expression matix:
> dim(exprs_mat)
[1] 112830   735
> exprs_mat[1:4, 1:3]
        Tarca_001_P1A01 Tarca_003_P1A03 Tarca_004_P1A04
1_at           6.062215        6.125023        5.875502
10_at          3.796484        3.805305        3.450245
100_at         5.849338        6.191562        6.550525
1000_at        3.567779        3.452524        3.316134

here is annotation data which contain experiment observation:
> dim(ano)
[1] 735   6
> head(ano)
                       SampleID   GA Batch     Set Train Platform
Tarca_001_P1A01 Tarca_001_P1A01 11.0     1 PRB_HTA     1    HTA20
Tarca_013_P1B01 Tarca_013_P1B01 15.3     1 PRB_HTA     1    HTA20
Tarca_025_P1C01 Tarca_025_P1C01 21.7     1 PRB_HTA     1    HTA20
Tarca_037_P1D01 Tarca_037_P1D01 26.7     1 PRB_HTA     1    HTA20
Tarca_049_P1E01 Tarca_049_P1E01 31.3     1 PRB_HTA     1    HTA20
Tarca_061_P1F01 Tarca_061_P1F01 32.1     1 PRB_HTA     1    HTA20

I intend to see how the genes in each sample are correlated with GA value of corresponding samples in the annotation data. How can I get my expected correlation matrix and filter out the genes by correlation value? any idea to make this happen correctly?
my attempt
library(limma)
fit <- limma::lmFit(exprs_mat, design = model.matrix( ~ 0 + t(ano$GA))
fit <- eBayes(fit)
topTable(fit, coef=2)

but this attempt gave me dimension error and I still didn't get highly correlated genes from exprs_mat with each ano$GAin anodata. I just removed 0+ and ran the code, then I got this error:
row dimension of design doesn't match column dimension of data object

I am pretty sure I did df[!is.na(df$GA), ], removed NA in both exprs_mat and ano metadata already, but the same error still persists. I also tried df[complete.cases(df),] to make sure I removed all NA values in my data.
How can I get those highly correlated genes correctly? any possible way to make this happen in R? any thought?
How can I cut down exprs_mat with the shortlist of genes? any idea?
desired output
I intend to filter out genes in the gene expression matrix expr_mat by using correlation value (only keep highly correlated genes), expected sub expression matrix should have same data structure as expr_mat. any way to get this done?

Have a look at the dimensions of what model.matrix() is returning, that should be telling.

The problem is that you are transposing the vector of GA values with t(ano$GA). Why would you do that? It produces a row matrix that is inappropriate for input to model.matrix.
All you need is
design <- model.matrix(~ano$GA)

and everything will work fine.


text:
creating population AF from set of individuals in 1000 genomes population?
What's the recommended way to create population-level Allele Frequency variant entries from a set of individuals (vcfs?) in the 1000 genomes (or similar) project?
For example, given the IBS individuals in the 1000 Genomes project:
http://www.internationalgenome.org/data-portal/population/IBS
I would like to create an IBS .vcf.gz that contains an IBS population-level AF entry for each of the variants present.
Any ideas on how to achieve this?
EDIT: this is related to my earlier question on how to create the most complete file for population allele frequencies, here:
what is the most complete vcf file of population allele frequencies that can be built/downloaded from public datasets?
I have seen that there is an Allele Frequency Calculator tool in EnsEMBL, which is at least available with GRCh37:
http://www.internationalgenome.org/allele-frequency-calculator-documentation
http://grch37.ensembl.org/Homo_sapiens/Tools/AlleleFrequency?db=core
There are also instructions on how to do this genome-wide using vcftools:
vcf-subset -c ~/1000g/IBS.samples.list ~/1000g/ALL.chr1.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz | fill-an-ac | bgzip -c > ~/1000g/IBS.chr1.20130502.vcf.gz

This seems to be generating correct files, and the AF can be calculated with a tool that combines AC and AN values, such as it is described here:
https://www.biostars.org/p/123468/
Any other ideas / warnings / pitfalls of the approach?

Are you looking to get this genome-wide? Or are you only interested in particular regions/variants?

Hi @Emily_Ensembl looking to do this genome-wide for as many populations as possible. I see there are 25-30 available in EnsEMBL from the 1000 genomes project.

bcftools has a plugin called fill-tags which should do the job for you easily, without having to mess around with lots of complicated bash scripting.
You could do something along the lines of this, presuming you have a .vcf file with all the IBS individuals in it:
bcftools +fill-tags in.vcf.gz -Oz -o out.vcf.gz

This will give each position in the .vcf file new INFO tags, one of which will be allele-frequency.
You can then pull off the frequency for a particular variant (CHR:POS) by:
bcftools query -t CHR:POS -f'%CHROM\t%POS\t%AF\n' out.vcf.gz

I will check it out. Thanks


text:
How to highlight specific cells in Seurat 2.4
I used Seurat 2.4 on our scRNA dataset to obtain the following tSNE plot.I was able to successfully extract cell IDs from the different clusters, and generate gene expression profiles. The analysis, and the biology makes sense. I have been however stuck in trying to highlight specific cells we are interested in using the Cell IDs (barcodes). For e.g, in this plot I am hoping to highlight specific cells (~ 120 of them) within Cluster 5 (blue cluster on the top). It can be either in featureplot mode or in this plot itself by an overlay, it doesn't matter. All I have to show are the 120 cells within the cluster.
For eg. if cluster 5 has cells atag-1 atgc-2 atat-3, cacc-4 cat-5... i want to recreate this plot using atag-1 atgc-2, cacc-4
I tried using do.identify = TRUE but that is a visual selection and can select regions, but is not useful to select specific cells. I was thinking if there is a way to subset and replace cluster 5 to include just the cells (~120 odd out of the total 200) but was unable to crack it. I tried by using cells.use but then realized that doesn't work either.Any ideas will be greatly appreciated, to get me out of this rut.
Thank you

Try to do this with PAGODA2 from Jean Fan: https://github.com/hms-dbmi/pagoda2

You can "highlight" any cells you'd like by creating a new column (or modifying an existing column) in the meta data slot within the Seurat object. From there you can modify the contents of that vector to your liking and plot the results.
It's not clear if you've been successful in obtaining the barcodes for the subpopulation to begin with, because that's a separate issue all together, and you need to clarify what attribute you are trying to subset the cells on because there are a number of ways to accomplish that. If you've already obtained the barcodes then I'll describe how you would plot the results below.
Let's assume you'd like to keep the original cluster IDs for all other cells. We'll start off by duplicating the cluster ID column (this isn't necessary, but it's good practice to preserve the original IDs).
n <- seurat_obj@meta.data # for convenience
n$new_column <- n$res.1.0

Now you need to modify the ids for cluster 5 to break them up into 2 populations.
n[rownames(n) %in% my_barcodes, "new_column"] <- "5b"

You can plot the new column using DimPlot with the argument group.by = "new_column", which will make the subpopulation stand out from the rest of the cells in cluster 5.


text:
"perl: warning: Setting locale failed." in RepeatMasker
I'm trying to run Repeatmasker in Linux on the command line with:
./RepeatMasker -species human -alu -gff -dir /mnt/lustre/users/Analysis/RepeatMaskerOutput /mnt/lustre/users/Software/STAR/bin/Linux_x86_64/HRG_Newest_and_Annotation/GCA_000001405.28_GRCh38.p13_genomic.fna

but I am getting this error:
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = "en_GB.utf8",
        LC_ALL = "",
        LC_CTYPE = "en_GB.utf8",
        LANG = "en_GB.utf8"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker version open-4.0.9
TRF program not configured! This version of RepeatMasker requires a
local installatoin of TRF.  Please visit:
               http://tandem.bu.edu/trf/trf.html
to obtain the current version.  Once installed please re-run the
RepeatMasker configure script to setup RepeatMasker to use the new
installation. at ./RepeatMasker line 442.

I have downloaded rmBlast and the latest TRF file. I have put the TRF file in my Users/username/bin PATH folder and used chmod a+x (and chmod 755) on it. I then changed my LD_LIBRARY_PATH using export to my glibc directory.
After this I used perl ./configure and directed the program to the PATH folder where the TRF file was. It gave a series of the same error messages as above but was 'successful' as shown here:

[user@nodec23() RepeatMasker]$ perl ./configure
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
RepeatMasker Configuration Program

**TRF PROGRAM**

  This is the full path to the TRF program.
 This is now used by RepeatMasker to mask simple repeats.

Enter path  [ /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/TRF.pm ]

Add a Search Engine:
   1. Crossmatch - Phrap             [ Configured ]
   2. RMBlast - Improved NCBI Blast  [ Configured, Default ]
   3. HMMER3.1 & DFAM                [ Un-configured ]
   4. ABBlast/WUBlast                [ Un-configured ]
   5. Done

Enter Selection: 5
 -- Setting perl interpreter...
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/Software/RepeatMasker/RepeatMasker/repeatmasker.help

I'm at a loss of what to do next. I have googled and tried changing each of the locale settings with export to the following:

en_GB
en_US.utf8
en_GB:2000
en_GB.utf8

but none have changed the error.
I am running a Windows 10 machine and using PuTTY to connect to an HPC which uses a Linux environment.
Locale settings are as follows:
locale: Cannot set LC_CTYPE to default locale: No such file or directory
locale: Cannot set LC_MESSAGES to default locale: No such file or directory
locale: Cannot set LC_ALL to default locale: No such file or directory
LANG=en_GB
LC_CTYPE="en_GB"
LC_NUMERIC="en_GB"
LC_TIME="en_GB"
LC_COLLATE="en_GB"
LC_MONETARY="en_GB"
LC_MESSAGES="en_GB"
LC_PAPER="en_GB"
LC_NAME="en_GB"
LC_ADDRESS="en_GB"
LC_TELEPHONE="en_GB"
LC_MEASUREMENT="en_GB"
LC_IDENTIFICATION="en_GB"
LC_ALL=

Also please note that I do not have administrator access.
Typing 'perl -v' gives the following:
perl: warning: Setting locale failed.
perl: warning: Please check that your locale settings:
        LANGUAGE = (unset),
        LC_ALL = (unset),
        LANG = "en_GB"
    are supported and installed on your system.
perl: warning: Falling back to the standard locale ("C").

This is perl, v5.10.1 (*) built for x86_64-linux-thread-multi

Copyright 1987-2009, Larry Wall

Perl may be copied only under the terms of either the Artistic License or the
GNU General Public License, which may be found in the Perl 5 source kit.

Complete documentation for Perl, including FAQ lists, should be found on
this system using "man perl" or "perldoc perl".  If you have access to the
Internet, point your browser at http://www.perl.org/, the Perl Home Page.

Typing "LC_ALL=C perl ./configure" and going throught the configurations leads to the same outcome but without the perl errors:
Enter Selection: 5
 -- Setting perl interpreter...

Congratulations!  RepeatMasker is now ready to use.
The program is installed with a the following repeat libraries:
  Dfam database version Dfam_3.0
  RepeatMasker Combined Database: Dfam_3.0
Further documentation on the program may be found here:
  /mnt/lustre/users/k1894998/Software/RepeatMasker/RepeatMasker/repeatmasker.help

However, the error remains when trying to run Repeatmasker.

I think it is quite simple you need to install perlbrew here
You then run the software under perlbrew it should work.
At present you are running the system perl5 and that is causing all sorts issues, for example the age of the installation and the fact system operations and your software are using the same version of perl.
All the instructions are in the website link.

Actually I'd recommend Anaconda/conda now.


text:
What is the formula for Mg values in TMM normalization for RNA Seq data?
I am reading through the paper "A scaling normalization method for differential expression analysis of RNA-seq data" by Mark D Robinson, Alicia Oshlack, available here.
In this paper they introduce a normalization method called "trimmed mean of M-values" (TMM), which is apparently a normalization used in the edgeR package.
What is confusing to me is that on page 3 they define the value
$M_g = \log_2\left(\frac{Y_{gk}/N_k}{Y_{gk'}/N_k'}\right)$
where $k$ is an experiment to be normalized and $k'$ is a reference experiment. A few pages down on page 7, however, they define
$M_{gk}^r = \frac{\log_2\left(Y_{gk}/N_k\right)}{\log_2\left(Y_{gr}/N_r'
\right)}$
where $k$ is an experiment to be normalized and $r$ is a reference experiment.
Are these two formulae intended to represent the same quantity? If they are supposed to be the same, which one has the correct order of logarithm/division? If they are not supposed to be the same, can anyone explain the difference? Thank you.