Ming Tang crazyhottommy

Three gotchas when using R for Genomic data analysis

During my daily work with R for genomic data analysis, I encountered several instances that R gives me some (bad) surprises.

some files such as bed file is 0 based. Two genomic regions:

	# How to test packages to be included in nixpkgs

	Thank you for helping us test packages to be included in nixpkgs!
	If you want to test a PR, first, you need to fork the nixpkgs
	Repository.

	To follow along, you first need to have Nix installed on your system.
	I recommend the installer from Determinate Systems, which you can find
	here: https://github.com/DeterminateSystems/nix-installer

	## add following to the end of ~/.bash_profile assuming ~/.bash_profile takes the precedence over ~/.profile (if present), ~/.bash_profile, or other bash configs.
	## for zsh or other shells, you may need to tweak following code.
	## Be careful "resetting" PATH variable and know the implications of PATH variable being either empty or missing critical system paths!

	################################### SET PATH ###################################
	## User specific environment and startup programs
	## Rewriting PATH ##
	## We want user and anaconda paths to precede system paths, esp. this string:
	## /usr/local/bin:/usr/bin:/usr/local/sbin:/usr/sbin which contains
	## most of system default libs,

	#!/bin/bash
	#
	# This script will find the consensus peak regions from peak files (in
	# BED format) of multiple samples by:
	#
	# 1. Converting the peak file of each sample into non-overlapping 3
	# cols BED file and concatenating them;
	#
	# 2. Sorting the concatenated file and Building a genome coverage
	# track in BedGraph, of which the value (the 3rd col) indicates the

	# code source for plotting
	# https://www.datanovia.com/en/blog/ggplot-title-subtitle-and-caption/
	# get source filename SO quest: http://stackoverflow.com/q/1815606/946850
	# get source fielname answer: https://stackoverflow.com/a/20282846
	# if you do use a hash, you will have to save these to a csv or write
	# a DB, perhaps an sqlite file that has the entry
	# for the source file and corresponding hash to keep track of the
	# hashes or just use filename and proj_dir,
	# and display that if you want to keep it simple,
	# open and reproducible

	--[[
	## Modulefile in lua syntax
	## Author: Samir Amin

	## Read about Lmod
	## https://lmod.readthedocs.io/en/latest/015_writing_modules.html
	## https://lmod.readthedocs.io/en/latest/050_lua_modulefiles.html
	## https://lmod.readthedocs.io/en/latest/020_advanced.html
	--]]

+(defun new-copyright ()
+  "Generate new copyright string."
+  (format-time-string "Copyright %Y Michael M. Hoffman <[email protected]>"))
+(define-auto-insert 'python-mode
+  `(
+    "Description: "
+    "#!/usr/bin/env python3.6" \n
+    "\"\"\"" (setq basename (buffer-file-basename)) ": " str \n
+    "\"\"\"" \n \n

	DATASETS = ["PBMC_8K", "PBMC_4k"]

	SALMON = "$BINS/salmon-0.14.0_linux_x86_64/bin/salmon"

	rule all:
	input: expand("quants/{dataset}/alevin/quants_mat.gz", dataset=DATASETS)

	rule salmon_quant:
	input:
	r1 = "reads/{sample}_1.fastq",

	# https://www.biostars.org/p/63225/

	FILE=VLI9_AA_S60_L008_R1_001.fastq.gz

	zcat $FILE \| head -n 40 \| awk '{if(NR%4==0) printf("%s",$0);}' \| od -A n -t u1 \| awk 'BEGIN{min=100;max=0;}{for(i=1;i<=NF;i++) {if($i>max) max=$i; if($i<min) min=$i;}}END{if(max<=74 && min<59) print "Phred+33"; else if(max>73 && min>=64) print "Phred+64"; else if(min>=59 && min<64 && max>73) print "Solexa+64"; else print "Unknown score encoding";}'