Genome Spot

CpG shores are the regions immediately flanking and up to 2 kbp away from CpG islands. These regions are interesting because methylation they are variably methylated in cancer and development (see reference ). Image courtesy of 7-themes.com To identify CpG shores, you first need a list of coordinates of CpG islands. If you have a linux computer with Emboss tools installed, you can determine CpG island positions with the cpgplot command . (if you don't have Emboss installed you can get CpG island coordinates from the UCSC table browser  just be wary of chromosome naming ie, "chr1" vs "1") cpgreport -score 17 -outfile genome.fa.cpg -outfeat /dev/stdout genome.fa | head Will produce a file like this one: ##gff-version 3 ##sequence-region 1 1 59999934 #!Date 2015-04-10 #!Type DNA #!Source-version EMBOSS 6.6.0.0 1 cpgreport sequence_feature 10469 11240 1299 + . ID=1.1 1 cpgreport sequence_feature 11284 11301 37 + . ID=1.2 1 cpgrepo...

When analyzing ChIP-seq data, we sometimes like to focus on promoters as a proxy of gene activation especially for histone marks that are normally associated with gene activation such as H3K9 acetylation or H3K4 tri-methylation. FeatureCounts has emerged as a competitor to HTSeq and BedTools MultiCov for counting reads across features (ie, exons, genes, promoters). FeatureCounts is great for RNA-seq because it can natively read GTF annotation files, but can't read BED format (that we use a lot in ChIP-seq analysis). In order to make featureCounts work, we need to extract the TSS coordinates and convert to a BED-like format that it can read ( SAF format ). In the below script, I extract the positions of the TSSs using a grep search followed by stripping only the neccessary information from the GTF, then using awk and BedTools to expand the region  around the TSS by 3kbp. #!/bin/bash #Generate an SAF file for TSS intervals from a GTF #Specify some parameters GTF=Homo_sapien...

Say you have a huge FASTA file such as genome build or cDNA library, how to you quickly extract just one or a few desired sequences? Use samtools faidx to extract a single FASTA entry  first index, then you can extract almost instantaneously. $ samtools faidx Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa real 0m37.422s $ time samtools faidx Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa MT real 0m0.003s Use bedtools getfasta extract portions of a FASTA entry Requires the regions of interest be in BED format. $ head Homo_sapiens.GRCh38_CpG_Islands.bed 1 10413 11207 1 28687 29634 1 51546 51882 1 121270 121549 The sequences of interest are extracted like this: $ bedtools getfasta -fi Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa -bed Homo_sapiens.GRCh38_CpG_Islands.bed -fo Homo_sapiens.GRCh38_CpG_Islands.fasta Make a blast database to access bunches of sequences quickly Note: you will need to download and install the BLAST+ package...

The unbiased nature of deep transcript sequencing makes it the ideal technology to discover novel uncharacterised genes. Lets screen our favourite RNA-seq experiment (azacitidine-treated AML3 human cells,  GSE55123 ) for novel expressed genes. I use Ensembl gene annotations. We'll start with preparing bed files of exons and gene bodies grep -w exon Homo_sapiens.GRCh38.76.gtf | tr '" ' '\t' \ | cut -f1,4,5,11 | uniq > Homo_sapiens.GRCh38.76_exons.bed grep -w gene Homo_sapiens.GRCh38.76.gtf | tr '" ' '\t' \ | cut -f1,4,5,11 > Homo_sapiens.GRCh38.76_genes.bed Now for convenience, I'll merge the data from the three replicates with samtools. samtools view -H SRR1171523_1.fastq.sort.bam > header.txt samtools merge -h header.txt Ctrl.bam SRR1171523_1.bam SRR1171523_2.bam SRR1171524_1.bam SRR1171524_2.bam SRR1171525_1.bam SRR1171525_2.bam samtools merge -h header.txt Aza.bam SRR1171526_1.bam SRR1171526_2.bam SRR1171527...

Counting the number of reads that align to certain genomic features is a key element of many next gen sequencing analysis pipelines. For RNA-seq, this is commonly used to count reads aligning to exons, while for ChIP-seq this is used to count reads over a promoter or other region of interest. There are several available tools for performing this task and in this post I will compare the three of the most commonly used: bedtools multicov htseq-count featureCounts I took one of the bam files from the recent RNA-seq series  of posts and subsampled it using samtools and shuf into file sizes of 1M, 2M, 5M, 10M reads, as well as the bam file containing 25.4M reads. I  then used the benchmarking script described in the previous post  to record execution time, CPU usage and peak memory for read counting to generate a gene-wise matrix. I used featureCounts in the single thread mode as well as the parallel mode (maximum of 8 cores). Execution time for read counting...

In this RNA-seq analysis series of posts , we're going through an RNA-seq analysis workflow. So far, we've downloaded, done QC and aligned the reads. In this post, we will count the reads aligning to exons and generate a count matrix for statistical analysis. My overall strategy is to use bedtools to count reads over exons. To do that, I'll need a bed file containing all exons. I could use UCSC genes, but I'd rather use the latest Ensembl version. Create exon bed file  To begin, we need a list of known exons. I will use the latest Ensembl human gene annotation gtf file from their ftp server. Then extract the exon entries with grep and then rearrange the gene accession and symbol information into a temporary bed file. The temporary bed file is exploded on gene name and then overlapping exons are merged with bedtools. #!/bin/bash wget ftp://ftp.ensembl.org/pub/release-76/gtf/homo_sapiens/Homo_sapiens.GRCh38.76.gtf.gz zcat Homo_sapiens.GRCh38.76.gtf.gz \ | gr...