6
Welcome to RNA-seq. This is not at all unusual, as RNA-seq raw data are not normally distributed. People typically approximate it via the Negative Binomial distribution. I suggest you check existing packages as Devon suggests in his comment for normalization and differential expression (or anything else). Typical choices are DESeq2, edgeR or limma, but there ...
6
Perhaps this small test script will help demonstrate some of the principles:
#!/usr/bin/env perl
use strict;
use warnings;
use Data::Dumper;
my @arr;
# define the current line
my $line = "foo\tbar\tbaz\n";
# ...
5
You're reinventing bedtools intersect (or bedops), for which there's already a convenient python module:
from pybedtools import BedTool
s3 = BedTool('s3.bed')
s4 = BedTool('s4.bed')
print(s4.intersect(s3, wa=True, wb=True, F=1))
The wb=True is equivalent to -wb with bedtools intersect on the command line. Similarly, F=1 is the same as -F 1.
3
This is a good task for Biopython:
from Bio import SeqIO
with open('input_file') as f:
motif = next(f).rstrip().lower()
for r in SeqIO.parse(f, 'fasta'):
print(r.id, ':', r.seq.lower().count_overlap(motif))
Note: in future I would recommend don't add the motif to the start of your FASTA file, keep it separate so as not to alter the ...
3
LAST has given the best results for me when I've tried to do this, although I agree with @user172818 that it's not a good idea to map really short reads. This is due to a combination of natural sequence duplication in long DNA sequences (e.g. see here), as well as abundant base calling differences present in single-molecule sequencing.
Minimising error is ...
3
Regarding your python code
If you want the experimental ranges that are entirely contained in one of the reference ranges, you need to have the coordinates in the following order:
cst <= tst < ten <= cen
If what you want are the experimental ranges that overlap one of the reference ranges, you need to have either the start or the end of the ...
3
Calculate in this context would just mean determine. You don't actually have to change every base to determine whether a change is silent. Instead, try building a hash table of amino acids to codons, which you can then quickly iterate over to find all silent changes.
BTW, try to avoid loaded terms like "mutation", even if your professor incorrectly used it. ...
3
In an ideal world, ribosomal RNA (as seen in your top hit) should be excluded from samples prior to sequencing. Where this is not possible (i.e. in the data that have been presented to you), it would be a good idea to exclude ribosomal genes (and any other common contaminants) prior to doing further analysis (including normalisation). I believe this is the ...
2
From what I can tell, the program output does not indicate whether the open reading frames (ORFs) are preceded by ribosome binding sites (RBSs). I'm not sure what parameters are used when establishing ORFs. To verify your RBS requirement, index the input DNA (using given ORF locations from biopython output, pay attention to strands, which are designated by 1/...
2
You should start by finding all instances of AAGGAGGTG and then look for ORFs that begin 6 to 9 nucleotides downstream of them. By begin, I mean find an AUG codon, and then extend until the next in-frame stop codon.
There's also no point in translating, just look for long stretches of sequence between AUG and UGA, UAA or UAG, those will define your reading ...
2
The following Perl documentation pages should be informative:
split - for splitting a scalar at all matches of a defined pattern
perlreftut - discusses the approach of anonymous variables and how to combine them
2
Let $n$ be the number of rows (in your case $n=5$) and $|O_j|$ the cardinality of $O_i$.
We can calculate the number of compatible columns as sum of the number of $O_i \subseteq O_j$, $O_j \subseteq O_i$ and $O_i, O_j$ disjoint.
number of $O_i \subseteq O_j$:
Here we need to find the number of supersets of our $O_i$. That means we need to find the number ...
2
What should I do with Ns (ambiguous bases), I've seen some people removing them or replacing them by gaps but this does not seem like the same information? Should they be left in?
As far as I understand, the lowercase letters stand for potential repeats. However, If I just ignore this, it is possible my sequences might not align to the repeat because it is ...
2
Suffix arrays vs Suffix tree is a problem quite often discussed in bioinformatics classes, but it's not that much used in practice. Most of people use already optimized mappers bowtie2 or bwa-mem.
I would advise you against your intention of mapping the reads using your own implementation of anything. You don't need to reinvent the wheel, just use an ...
2
Are you not allowed to use the web BLAST tool? That's what I use if I need to quickly find out the likely origin of a DNA sequence.
There are command-line ways to do the same thing, but I only use them if I need to search hundreds or thousands of sequences at once.
2
Unless it's not obvious, the problem has no unique solution. An obvious example is a string which uses entirely one character; the problem is then to determine its length, which you can only do statistically.
As you've correctly surmised, this is more or less the short-read de novo assembly problem, with one modification: the "read errors" are all ...
2
I like using DESeq2. There's a great document written by the developers about how to process gene count tables and look at differential expression:
https://bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html
The document also mentions other tools that have similar approaches:
Other Bioconductor packages with similar aims are edgeR, ...
1
[Comments from other post migrated as answer]
The parenthetical statement in the first bullet point says without corrections. In which case, if an unseen nucleotide has a probability of zero, the probability of the sequence is zero. If the unseen nucleotides had the probability of NaN, then the sequence has a probability of NaN, which is a more thorough ...
1
As mentioned by Greg, list of short-read sequencing platforms is listed on wikipedia.
Of those you mentioned Illumina and Ion Torrent are short-read sequencing platforms.
Usually what makes short reads short is a PCR step in the library prep or directly on the flowcell to amplify the signal, which is the limitation for how long stretches of DNA are possible ...
1
You will want to create a self loop for a $k$-mer $s$ whenever $s_1,...,s_{k-1} = s_2,...,s_k$. This will only happen when all of the characters of the $k$-mer are identical, as $s_1=s_2, s_2=s_3,...,s_{k-1}=s_k$.
1
Without giving away too much for your homework check out the sample() function in R. Think about how you can pass DNA bases to this function to generate random sequences.
1
Thanks for all your answers.
Although the pipeline packages would have been useful I was trying to normalise the data before analysing it later and so I used the bestnormalize package which looked like it did a pretty good job.
1
Another approach, that does not replace the methods already mentions, but could be useful if you don't have raw counts to start with, or if you want to screen hundreds of genes, and you want to focus on the most interesting.
either log transform or zscale your expression values
do a boxplot to see how your gene is expressed in the different tissues
run a ...
1
Your alignment is
SDRVIKAAIFDPIQPDF---G-----PVYFGLGHVH
RDLVERLFILDMI-PGLIKAGDSFPIPVALMINHIF
For each position in the alignment you calculate the score for that alignment. For position 1 we'd look up S vs R in the matrix and find a score of -1. You continue doing this until you hit the first -, which is not in the matrix. The first - is a gapopening, each ...
1
According to the book you mentionned Introduction to Genomics SecondEdition by Arthur MLesk p.113
(a) What fraction of a genome could you expect to assemble from eightfold coverage?
(b) What total gap length would you expect in an assembly of a 2 Mb target genome size from eightfold coverage?
(c) How many gaps would you expect in an assembly of a 2 Mb ...
1
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 ...
1
I tried to fix the indentation of the code you posted, a typo in some variable names, plus a few other issues that caused it to not work as expected (see the comment in the code):
#!/usr/bin/env python3
import sys
sequence = sys.stdin.readlines()
d = {}
temp_genename = None
temp_sequence = None
# The lines returned by readlines() contain end-of-line ...
1
A protein domain is a conserved part of a given protein sequence and (tertiary) structure that can evolve, function, and exist independently of the rest of the protein chain.
This is different from global protein similarity, so similarity to neurotoxin type F doesn't necessarily mean that sequence A has similar domain to your known toxin gene. Your blast ...
1
pyranges answer:
# pip install pyranges
# or
# conda install -c bioconda pyranges
import pyranges as pr
s3 = pr.read_bed("s3.bed")
s4 = pr.read_bed("s4.bed")
s3.intersect(s4, how="containment")
Answer with setup (Edit):
import pandas as pd
from io import StringIO
import pyranges as pr
c1 = """1 10 20
1 5 20
2 20 30
2 25 30
1 10 50
2 ...
1
You could use BEDOPS, instead:
$ sort-bed s3.txt > s3.bed
$ sort-bed s4.txt > s4.bed
$ bedops --element-of 1 s3.bed s4.bed > answer.bed
If you need to run it from within Python, you could use subprocess.check_output():
import subprocess
...
try:
result = subprocess.check_output("bedops --element-of 1 %s %s > %s" % (set_a_fn, set_b_fn, ...
Only top voted, non community-wiki answers of a minimum length are eligible