This Nextflow[1] based workflow allows alignment of either single or paired ends reads to a reference transcriptome described as a genome in fasta file plus an annotation in GTF format using STAR[2] aligner considering the variant information provided as a VCF file.
STAR implements the WASP[3] method for filtering of allele specific alignments and reports within the resulting alignments which variant is found and from which allele. A python script based on pysam [4] efficiently separates the alignment generated by STAR into four different files:
- Allele A
- Allele B
- Neither A nor B
- Ambiguous (found both variants on one or different pairs)
HTseq-count[5] tool is then used for counting tags separated in those three categories mapping to the genes. Finally a report is generated using multiQC[6] that summarize the results of each step together with an initial QC evaluation fo raw reads done using FastQC[7]
You need to install either Docker or Singularity and Nextflow.
Then you can clone the repository:
git clone --depth 1 [email protected]:biocorecrg/allele_specific_RNAseq.gitor
git clone --depth 1 https://github.com/biocorecrg/allele_specific_RNAseq.git
depending on your GitHub configuration.
You need to extract the SNP information from the global VCF file, so first of all you need to download the VCF file with the index:
wget ftp://ftp-mouse.sanger.ac.uk/REL-1505-SNPs_Indels/mgp.v5.merged.snps_all.dbSNP142.vcf.gz
wget ftp://ftp-mouse.sanger.ac.uk/REL-1505-SNPs_Indels/mgp.v5.merged.snps_all.dbSNP142.vcf.gz.tbithen the reference genome:
wget ftp://ftp-mouse.sanger.ac.uk/ref/GRCm38_68.fa
and the annotation from Ensembl. We used the version Mus_musculus.GRCm38.68 not available in Ensembl archive.
You can use the pipelines specifying between Docker and Singularity containers by using the options
-with-singularity
or
-with-docker
The module makeAnno can be used for generating a VCF file with SNP for the interesting species and a genome with SNP position masked with Ns.
For using the module:
cd makeAnno
NXF_VER=20.01.0 nextflow run make_anno.nf -with-singularity -bg --vcffile mgp.v5.merged.snps_all.dbSNP142.vcf.gz --speciesA CAST_EiJ --speciesB 129S1_SvImJ --genome GRCm38_68.fa --outvcf CAST_EiJ-129S1_SvImJ.vcf > logOr in case you want to compare against the reference genome:
cd makeAnno
NXF_VER=20.01.0 nextflow run make_anno.nf -with-singularity -bg --vcffile mgp.v5.merged.snps_all.dbSNP142.vcf.gz --speciesA CAST_EiJ --speciesB reference --genome GRCm38_68.fa --outvcf CAST_EiJ-C57BL_6J.vcf > logThis can take some memory for masking the genome. In the output folder named filteredVCF you will find both gzipped masked genome and gzipped vcf annotations.
cd allele_specific_RNAseq;
NXF_VER=20.01.0 nextflow run as_rnaseq.nf -with-singularity -bg > logoptionally you might want to check your pipeline on Nextflow's Tower website. You need to register using an istitutional mail and set the token provided in a variable as described:
export TOWER_ACCESS_TOKEN=<<<<<TOKEN NUMBER>>>>>>
NXF_VER=20.01.0 nextflow run as_rnaseq.nf -with-singularity -bg -with-tower > log you can check the status of your pipeline live on the tower website.
The parameters for running the pipeline are defined in the file params.config that can be changed accordingly.
| parameter | value |
|---|---|
| reads | $baseDir/test/*_{1,2}.fastq.gz |
| genome | $baseDir/test/GRCm38_68_19.masked.fa.gz |
| annotation | $baseDir/test/Mus_musculus.GRCm38.68_19.gtf |
| strandness | reverse |
| variants | $baseDir/test/19_filt.vcf.gz |
| output | $baseDir/output_test |
| single | NO |
| varcut | 1 |
| title | Allele specific RNAseq project |
| subtitle | This is my wonderful RNA experiment |
| PI | Luca Cozzuto |
| User | Luca Cozzuto |
| UCSCgenomeID | mm10 |
| [email protected] |
providing a real email address will deliver a mail with the multiqc report when the analysis is finished.
Fastq paired ends reads can be either plain or gzipped. You need to specify this syntax for capturing the identifier from single end samples:
PATH/*.fastq.gz
and this for paired end ones:
PATH/*_{1,2}.fastq.gz
Gzipped masked fasta file of the genome obtained running makeAnno
It can be either "unstranded" or "forward" or "reverse". It depends on the kind of sequencing.
GTF file
Gzipped VCF file obtained running makeAnno
YES: single end reads. NO: paired ends
Number of SNP needed for assigning a read to a variant.
The following folder will contain the final outputs:
- Index: the indexed genome
- QC: containing fastQC results
- Alignments: containing sorted BAM files as they were from normal bulk RNAseq
- Report: A detailed report for the pipeline run, that will be sent via email.
- Allele_alignments: containing BAM files containing reads with variants marked as: alleleA, alleleB, ref and ambiguous. The count is done considering the strand protocol used and specified as a parameter.
- cut_N * A folder that is generated for each SNP cut off chosen containing the following sub-folders:
- Allele_single_counts: containing the count per gene per sample and single alleles. The count is done considering the strand protocol used and specified as a parameter.
- Allele_merged_Counts: containing the count per gene per allele in a single file for each sample. The count is done considering the strand protocol used and specified as a parameter.
- Counts: composite counts per gene per sample. No distincion between alleles. Four fields: gene id, tot counts considering the sequencing unstranded, tot counts considering the sequencind done on the forward strand, tot counts considering reverse strand. This information is useful for validating the strand specific protocol used in the sequencing step.
- Proportions: for each sample the read count per allele are divided for the sum of counts of the variants and multiplied for the composite count. The count is done considering the strand protocol used and specified as a parameter.
propA = alleleA/(alleleA+alleleB) * composite
propB = alleleB/(alleleA+alleleB) * composite
React
Typescript
Django
Laravel
Facebook
Microsoft
Google
Alibaba
D3
Tencent