<jupyter_start><jupyter_text>Start up an AWS ubuntu instance and install docker. For this analysis, the following specs should be more than enough to run the bioinformatics pipeline.<jupyter_code>docker build  .
docker push<jupyter_output><empty_output><jupyter_text>Pull the docker image.<jupyter_code>sudo docker pull jonessarae/nlm_workshop:seq_tools<jupyter_output><empty_output><jupyter_text>Download SRA files with sra_download.sh (check permissions) Already hardcoded with the SRA accession number.<jupyter_code>./sra_download.sh<jupyter_output><empty_output><jupyter_text>Run the docker container and add home directory to container. The rest of this notebook will be done within the container.---<jupyter_code>sudo docker run -it -v ~:/mnt jonessarae/nlm_workshop:seq_tools<jupyter_output><empty_output><jupyter_text>Convert SRA files into FASTQ files and split the files into read 1 and read 2. The command fastq-dump for SRA-toolkit does not use multi-threading. Better to find another package that can multi-thread. Can avoid splitting? In BWA can use -p option. How about others?<jupyter_code>fastq-dump --split-files --origfmt --gzip SRR2989954.sra
gunzip SRR2989954_1.fastq.gz SRR2989954_1.fastq.gz<jupyter_output><empty_output><jupyter_text>Download reference genome hg19.2bit and tool twoBitToFa, convert 2bit file to fasta file, and index it with bwa. Use -t to add more cores for running bwa index faster.<jupyter_code>wget http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/hg19.2bit

rsync -aP rsync://hgdownload.soe.ucsc.edu/genome/admin/exe/linux.x86_64/twoBitToFa
chmod 744 twoBitToFa
./twoBitToFa hg19.2bit hg19.fa

bwa index -a bwtsw hg19.fa<jupyter_output><empty_output><jupyter_text>Download reference genome hg19 annotation file (GTF). There are two files to choose from. Not sure which to use. https://www.gencodegenes.org/releases/19.html<jupyter_code>wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz<jupyter_output><empty_output><jupyter_text>Exome: Align paired reads to the reference genome and convert it to a bam file. With an AWS instance with 96 cores, this step took ~25 minutes (on file 54) with 96 threads.<jupyter_code>bwa mem -t 8 ref/hg19.fa SRR2989954_1.fastq SRR2989954_2.fastq | samtools view -Sb - > SRR2989954.bam<jupyter_output><empty_output><jupyter_text>Sort the bam files by genomic order.<jupyter_code>samtools sort -@ 4 SRR2989954.bam -o SRR2989954_sorted.bam
samtools sort -@ 8 SRR2989963.bam -o SRR2989963_sorted.bam<jupyter_output><empty_output><jupyter_text>RNA: Make genome indices that STAR requires. Always ensure read/write priveleges on files/folders. Check what read length should be: ReadLength-1. Does it matter if not 99? Change number of threads to 8. Mention spec requirements.<jupyter_code>mkdir star_indices
STAR --runThreadN 4 --runMode genomeGenerate --genomeDir star_indices --genomeFastaFiles ref/hg19.fa --sjdbGTFfile gtf//gencode.v19.annotation.gtf --sjdbOverhang 100<jupyter_output><empty_output><jupyter_text>RNA: Align paired reads to the reference genome. Don't forget to make a mapped_rna directory.<jupyter_code>mkdir mapped_rna
STAR --readFilesIn SRR2989969_1.fastq SRR2989969_2.fastq --runThreadN 8 --genomeDir star_indices --outFileNamePrefix mapped_rna/SRR29899692pass --genomeLoad NoSharedMemory --sjdbGTFfile gtf/gencode.v19.annotation.gtf --outSAMtype BAM SortedByCoordinate --twopassMode Basic<jupyter_output><empty_output><jupyter_text>Exome: GATK4 from bioconda.<jupyter_code>java -Xmx4g -jar GenomeAnalysisTK.jar -I SRR2989954_sorted.bam -R ref/hg19.fa -T RealignerTargetCreator -o SRR2989954.intervals –known Mills_and_1000G_gold_standard.indels.hg19.vcf --read_filter MappingQualityZero 
java -Xmx4g -jar GenomeAnalysisTK.jar -I SRR2989954_sorted.bam -R ref/hg19.fa -T IndelRealigner -targetIntervals sample.intervals -o SRR2989954_realign.bam -known Mills_and_1000G_gold_standard.indels.hg19.vcf --read_filter MappingQualityZero

java -Xmx4g -jar GenomeAnalysisTK.jar -I SRR2989963_sorted.bam -R ref/hg19.fa -T RealignerTargetCreator -o SRR2989954.intervals –known Mills_and_1000G_gold_standard.indels.hg19.vcf --read_filter MappingQualityZero 
java -Xmx4g -jar GenomeAnalysisTK.jar -I SRR2989963_sorted.bam -R ref/hg19.fa -T IndelRealigner -targetIntervals SRR2989963.intervals -o SRR2989963_realign.bam -known Mills_and_1000G_gold_standard.indels.hg19.vcf --read_filter MappingQualityZero<jupyter_output><empty_output><jupyter_text>Filter out reads that have mapping quality of < 20.<jupyter_code>samtools view -b -q 20 SRR2989954_realign.bam > SRR2989954_m20.bam
samtools view -b -q 20 SRR2989963_realign.bam > SRR2989963_m20.bam
samtools view -b -q 20 mapped_rna/SRR29899692passAligned.sortedByCoord.out.bam > SRR2989969_m20.bam<jupyter_output><empty_output><jupyter_text>Deduplicate files. We installed picard separately, but picard is also part of GATK4.<jupyter_code>java -jar picard.jar MarkDuplicates I=SRR2989954_m20.bam O=SRR2989954_m20dedup.bam REMOVE_DUPLICATES=true METRICS_FILE=metrics.txt
java -jar picard.jar MarkDuplicates I=SRR2989963_m20.bam O=SRR2989963_m20dedup.bam REMOVE_DUPLICATES=true METRICS_FILE=metrics.txt
picard MarkDuplicates I=SRR2989969_m20.bam O=SRR2989969_m20dedup.bam REMOVE_DUPLICATES=true METRICS_FILE=metrics.txt<jupyter_output><empty_output><jupyter_text>Generate VCF for multiple BAM files. Not sure how they run this except that it outputs vcf format. Not sure where the normal sample is in the VCF file. Assumed they combined all three with normal.<jupyter_code>samtools mpileup -f hg.19a SRR2989954_m20dedup.bam SRR2989963_m20dedup.bam SRR2989969_m20dedup.bam -v -o p1_ov.vcf<jupyter_output><empty_output><jupyter_text>Calls variants from a mpileup dataset and produces a VCF.<jupyter_code>varscan mpileup2snp p1_ov.vcf --min-coverage 5 --min-reads2 0 --min-avg-qual 20 --min-var-freq 0 --output-vcf 13<jupyter_output><empty_output><jupyter_text>Run snpEff, a variant annotation and effect prediction tool. It annotates and predicts the effects of variants on genes. Not sure which version of the GRCh37 database was used.<jupyter_code>#snpEff databases | grep -i sapiens
snpEff download GRCh37.75
snpEff GRCh37.75 p1_ov.vcf > p1_ov.ann.vcf<jupyter_output><empty_output><jupyter_text>Count total reads and compare to supplementary figure 1.<jupyter_code>samtools flagstat <sample.bam><jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.08044901777362021, 'fraction_numerical': 0.0661053944496414, 'mean_word_length': 6.096238938053097, 'pattern_counts': {'":': 0, '<': 75, '<?xml version=': 0, '>': 79, 'https://': 1, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '18', 'dir': 'jupyter-structured-clean-dedup', 'id': '164810', 'license': 'MIT', 'n_tokens_mistral': 2665, 'n_tokens_neox': 2323, 'n_words': 573, 'path': 'nlm_workshop_ovarian_cancer.ipynb', 'repo_name': 'NLM-Reproducibility-Project/Repro_PMID_26735499-', 'size': '6414'}