This is an RNA editing pipeline that calls RNA editing events from RNA-seq data. It uses the human genome as reference, but can easily be modified to use any reference genome. The aim of this repository is to provide the code used in our corresponding article.
-
Clone the repository:
git clone https://github.com/oscar-franzen/rnaed/
-
Install software dependencies. The tools should be in your
$PATH
.
- GNU Wget
- STAR v. 2.5.1b
- samtools v. 1.3
- java v. 1.8.0
- Picard v. 1.112
- bamtools
- bedtools
- GNU Awk
- BioPerl (Bio::SeqIO)
- Perl 5.x
- gmap-2016-05-25
- NCBI liftOver
- GNU parallel
Optional:
This bash script reports missing dependencies (except Picard):
chmod +x check_missing.sh
./check_missing.sh
-
Prepare a file called "ucsc_146.snps.GRCh38.final", goto http://genome.ucsc.edu/cgi-bin/hgTables and select group "Variation" and track "All SNPs(146)" and Assembly "Dec. 2013(GRCh38/hg38)". In output format, choose "BED". Type in the filename ucsc_146.bed and press "Get output".
Then execute:
awk -F '\t' '{print $1"\t"$3}' ucsc_146.bed > ucsc_146.snps.GRCh38.final
-
Repeat step 2, but instead for version 141. The final file should be called ucsc_141.snps.GRCh38.final.
-
Register and download COSMIC. The license of this database is not compatible with unrestricted distribution. Only the VCF files are needed. Put all files in a directory called
COSMIC
. Execute:for file in ./COSMIC/Cosmic* do cat ${file} | grep -v -E '^#' | awk -F '\t' '{print $1"\t"$2}' > ${file}.C done
-
Prepare a file named
mapping_bias.final.txt
. This file should have a one line header and four space-separated columns:
- the sample name (taken from the fastq file:
sample_name.fastq
) - number of reads mapping mapping in the forward direction on protein-coding exons
- number of reads mapping in the reverse direction on protein-coding exons
- sequence read length
One way to retrieve the counts is to use the program htseq. 7. Put your *.fastq files in a directory called `fastq'. 8. Run:
chmod +x pipeline.sh *.pl
./pipeline.sh
- If everything worked, the called event will be in
./calls/
. One file for each fastq file. Each output file has has the following columns:
- chromosome
- position
- editing depth
- reference base
- edited base
- the pileup
- total depth
- position in sequencing read. The format is position:number of reads.
- strand. Only relevant for strand-specific sequencing.
- overlapping annotation(s) (e.g., gene, repeat)
- Franzén O, Ermel R, Sukhavasi K, Jain R, Jain A, Betsholtz C, Giannarelli C, Kovacic JC, Ruusalepp A, Skogsberg J, Hao K, Schadt EE, Björkegren JLM. (2018) Global analysis of A-to-I RNA editing reveals association with common disease variants. PeerJ 6:e4466 https://doi.org/10.7717/peerj.4466
- Oscar Franzén, [email protected]