Taking a different tack from other answers, there's lots of tools for pipelines in Python. Note: there was a time when people would use "pipeline" to refer to a shell script. I'm talking about something more sophisticated that helps you decompose an analysis into parts and runs it robustly.
Snakemake is my favourite. It's (nearly) pure Python and can ...
An epic question. Unfortunately, the short answer is: no, there are no widely used solutions.
For several thousand samples, BCF2, the binary representation of VCF, should work well. I don't see the need of new tools at this scale. For a larger sample size, ExAC people are using spark-based hail. It keeps all per-sample annotations (like GL, GQ and DP) in ...
The main difficulty here is the use of GRCh38. Unfortunately, despite the fact that it's more than four years old, the major ethnicity-labeled public datasets (1000 Genomes, gnomAD when allele frequencies are enough) still aren't available for that reference. It is necessary to perform a liftover operation, or just use overlapping rsIDs and hope for the ...
According to the SAM specification, the 3rd field of a SAM line (RNAME) is:
RNAME: Reference sequence NAME of the alignment. If @SQ header lines
are present, RNAME (if not ‘*’) must be present in one of the SQ-SN
tag. An unmapped segment without coordinate has a ‘*’ at this field.
However, an unmapped segment may also have an ordinary coordinate ...
The FAQ offers an answer:
I'm setting up my own Blat server and would like to use the same parameter values that the UCSC web-based Blat server uses.
We almost always expect there to be some small differences between the hgBlat/gfServer and the stand-alone command-line blat. The best matches can be found using pslReps and pslCDnaFilter utilities. ...
I've been using the LoFreq* caller for exactly this. It is designed to find variants with very low frequency, so is well suited for this type of analysis.
LoFreq* (i.e. LoFreq version 2) is a fast and sensitive variant-caller for inferring SNVs and indels from next-generation sequencing data. It makes full use of base-call qualities and other sources of ...
There is a recent paper that attempts to do this:
ISOWN: accurate somatic mutation identification in the absence of
normal tissue controls.
In this work, we describe the development, implementation, and
validation of ISOWN, an accurate algorithm for predicting somatic
mutations in cancer tissues in the absence of matching normal tissues.
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.
BioPython has some good tools for processing reads and alignments.
There is a python library wrapping samtools so many of the samtools calls can be used directly as python objects and calls
I would use subprocess to call the aligner and specify the output to a bam ...
Just use bcftools view for filtering:
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz > output.vcf
To truncate this list to 65,000 SNPs count the header lines, sum them to the number of SNPs you like to have and use head -n.
$ bcftools view -h input.vcf.gz|wc -l
$ bcftools view -i 'AF>0.3 && AF<0.7' input.vcf.gz | ...
What is "S13"?
S# is most likely the sample number. If so, you do not want to combine these.
What is "L002_l1"?
L00# is the lane number (of the 8 Illumina flow cell lanes). 'I1_' (not 'L') are index files and are not needed (see also https://www.biostars.org/p/374581/).
What is "L002_R1"?
R1 is one end of paired-end sequencing, R2 is the other end. ...
If you're just doing alignment and conversion to a sorted BAM file, there's no need to run it through python. A simple pipe on the unix command line works just as well (and probably runs faster):
hisat2 -x genome.index -1 reads_R1.fq.gz -2 reads_R2.fq.gz |
samtools sort > reads_vs_genome.bam
What's your reference for the definition of STRs? I think it is still ambiguous among the community. Wikipedia states motifs of >=2 bp. However other references include homopolymers as well:
Source 1, source 2 and source 3.
I agree that using a specialized tool is probably a good idea.
Nevertheless if you want to stick with Python, I suggest using plumbum instead of subprocess. It has a very nice syntax for this kind of problems.
from plumbum import local
# load a command
samtools = local['samtools']
# you can then easily redirect
(samtools['view', '-bS', '../some.sam'] >...
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. ...
Here is one definition of the classical MHC region: https://www.ncbi.nlm.nih.gov/grc/human/regions/MHC?asm=GRCh38.p13
which defines it as the 5mb region chr6:28510120-33480577 in GRCh38 coordinates.
The "extended MHC" defined here  is 7.2 mb and encompasses an additional "extended Class I" region upstream of the region given above. The ...
It turns out there is a tool that does this beautifully. I used the proteinToGenome function from the ensembldb package.
The code looks like this:
# Define the Ensemble database object for the
edbX <- filter(EnsDb.Hsapiens.v75, filter = ~ seq_name == "X")
# Set the first and last amino acid you want to convert and the ensembl protein id
tldr: Remove the $ from the command.
I imagine you're literally typing $hisat2, where you mean to instead type hisat2. The $ is just meant to show the end of the command prompt. If you actually type $hisat2, you're not executing the hisat2 program, but whatever the $hisat2 variable is set to. Since it's likely not set to anything, it's just ignored. For ...
The Bioinformatics Stronghold at Rosalind is an excellent resource for building your understanding of many concepts and algorithms you'll see frequently in bioinformatics and genomics.
I'd also say that, if you can make it work, an internship where you can "get your hands dirty" with real problems is a very valuable experience. This option may not be ...
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 ...
I finally found another answer to my question. Please read this great article in May 12 2017 BioMed Central (BMC) Bioinformatics article titled Ranking metrics in gene set enrichment analysis: do they matter?.
Also, please read this blog , Diving into Genetics and Genomics: Gene Set Enrichment Analysis (GSEA) explained.
After reading these two articles, my ...
Using bcftools is the right way, but you could also greatly simplify the command you're using. You can do the same thing with a single zgrep call:
zgrep -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz
And to limit to the first 65000 results:
zgrep -m 65000 -P '^[^#]+\t\d+\t\S+\t\w\t\w\t.*;AF=0\.[3-7]' my.vcf.gz
Finally, to keep the header of the ...
It seems so, yes.
UCSC uses a 10Kb definition and representation for telomeres in human. In the table you shared this is consistent across all chromosomes, not only chr21. This 10kb representation is an approximation to the mean value observed in newborns and is supported by evidence like the one linked below.
Taken from Wikipedia:
A telomere is a ...
Try with PsyGeNet:
PsyGeNET (Psychiatric disorders Gene association NETwork) is a resource for the exploratory analysis of psychiatric diseases and their associated genes
It is based on literature mining and curated by experts.
Here is the publication if you want more details:
PsyGeNET: a knowledge platform on psychiatric disorders and their genes
Alright, so there are a number of problematic patterns in your code - as far as I understand what you are trying to do. Next time, try to post a reproducible example that people can use and more people will be willing to help.
combination_labels = 
combination_counts = 
for lineage in lineages:
Declaring these two lists before the loop, then ...