RatGTEx pipeline

STAR --runMode genomeGenerate \
    --genomeDir {params.outdir} \
    --genomeFastaFiles {input.fasta} \
    --sjdbGTFfile {input.gtf} \
    --sjdbOverhang {params.overhang} \
    --runThreadN {resources.cpus}

bcftools view -s {wildcards.rat_id} --min-ac=1 -O z -o {output} {input}

STAR --runMode alignReads \
        --runThreadN {resources.cpus} \
        --genomeDir {params.index_dir} \
        --readFilesIn {params.fastq_list} \
        --readFilesCommand zcat \
        --twopassMode Basic \
        --varVCFfile <(zcat {input.vcf}) \
        --waspOutputMode SAMtag \
        --quantMode TranscriptomeSAM \
        --outSAMstrandField intronMotif \
        --outSAMattrRGline {params.read_groups} \
        --outSAMtype BAM SortedByCoordinate \
        --outSAMunmapped Within \
        --outFileNamePrefix {params.prefix}

grep -v '^#' {input} | awk '$3=="exon"' | cut -f1,4,5 | gzip -c > {output}

bcftools view {input.vcf} \
    --regions-file {input.regions} \
    -Ou | bcftools view \
    --min-af 0.2:minor \
    -i 'F_MISSING<0.1' \
    -Oz -o {output.vcf}

tabix -p vcf {output.vcf}

gatk ASEReadCounter \
    -R {input.ref} \
    -I {input.bam} \
    -V {input.vcf} \
    -O {output} \
    --min-depth-of-non-filtered-base 10 \
    --min-mapping-quality 250 \
    --min-base-quality 15

rsem-prepare-reference \
    {input.fasta} \
    ref/rsem_reference/rsem_reference \
    --gtf {input.gtf} \
    --num-threads {resources.cpus}

rsem-calculate-expression \
    {params.paired_end_flag} \
    --num-threads {resources.cpus} \
    --quiet \
    --estimate-rspd \
    --no-bam-output \
    --alignments {input.bam} \
    {params.ref_prefix} \
    {params.out_prefix}

gzip {params.out_prefix}.genes.results

python3 src/collapse_annotation.py \
    ref/Rattus_norvegicus.Rnor_6.0.99.gtf \
    ref/Rattus_norvegicus.Rnor_6.0.99.genes.gtf

python3 src/assemble_expression.py {params.rsem_dir} {input.anno} {params.prefix}

bgzip {params.prefix}.log2.bed
bgzip {params.prefix}.tpm.bed
bgzip {params.prefix}.iqn.bed
bgzip {params.prefix}.iqn.filtered.bed

tabix {params.prefix}.log2.bed.gz
tabix {params.prefix}.tpm.bed.gz
tabix {params.prefix}.iqn.bed.gz
tabix {params.prefix}.iqn.filtered.bed.gz

plink2 --make-bed \
    --vcf {input.vcf} \
    --keep {input.samples} \
    --maf 0.01 \
    --mac 2 \
    --max-alleles 2 \
    --out {params.prefix}

plink2 \
    --bfile {params.prefix} \
    --geno 0.05 \
    --maf 0.05 \
    --indep-pairwise 200 100 0.1 \
    --out {params.pruned_prefix}

plink2 \
    --bfile {params.prefix} \
    --extract {params.pruned_prefix}.prune.in \
    --export vcf bgz id-paste=iid \
    --out {params.pruned_prefix}

Rscript src/covariates.R {input.vcf} {input.bed} {params.n_geno_pcs} {params.n_expr_pcs} {output}

python3 src/run_tensorqtl.py \
    {params.geno_prefix} \
    {input.bed} \
    {output} \
    --covariates {input.covar} \
    --mode cis

python3 src/run_tensorqtl.py \
    {params.geno_prefix} \
    {input.bed} \
    {output} \
    --covariates {input.covar} \
    --cis_output {input.cis} \
    --mode cis_independent

python3 -m tensorqtl \
    {params.geno_prefix} \
    {input.bed} \
    {wildcards.tissue} \
    --covariates {input.covar} \
    --output_dir {params.outdir} \
    --mode cis_nominal

python3 -m tensorqtl \
    {params.geno_prefix} \
    {input.bed} \
    {wildcards.tissue} \
    --covariates {input.covar} \
    --output_dir {params.outdir} \
    --output_text \
    --batch_size 10000 \
    --mode trans

python3 src/tensorqtl_all_signif.py {input.perm} {params.nom_prefix} {output}

python3 src/tensorqtl_all_cis_pvals.py {params.nom_dir} {output}

python3 tools/aFC/aFC.py \
    --vcf {input.vcf} \
    --pheno {input.bed} \
    --qtl <(python3 src/prepare_qtl_for_afc.py {input.qtl} {input.qtl_indep}) \
    --cov {input.covar} \
    --log_xform 1 \
    --output {output}