For example, to create a paired end same file with fully formatted read groups, you can run:
bwa sampe -i MARK -m DEPRISTO -l LIBRARY -p ILLUMINA -h Homo_sapiens_assembly18.fasta first.fastq.sai secondfastq.sai first.fastq second.fastq > pe.sam
This will result in a file that looks like:
@HD VN:1.0 SO:unsorted
@SQ SN:chrM LN:16571
@SQ SN:chr1 LN:247249719
@SQ SN:chr2 LN:242951149
[More @SQ lines removed for clarity]
@RG ID:MARK SM:DEPRISTO LB:LIBRARY PL:ILLUMINA
61CC3AAXX100125:6:48:11631:21212 81 chrM 8 37 76M chr17 21944837 0 GGTCTATCACCCTATTAACCACTCACGGGAGCTCTCCATGCATTTGGTATTTTCGTCTGGGGGGTATGCACGCGAT 99:<->=D>IF7.:7,H
SNP and MicroIndel detection
Novoalign detects mismatches and short indels in short reads from next-generation sequencing platforms. The recipes below contain command line usage to generate list of single nucleotide polymorphisms SNPs and Indels for further analysis.
Generate sorted alignments in BAM format
Inititally we need to align our file(s) to the reference genome:
novoalign -d hg18.nix -f reads1.fastq -o SAM 2> stats.txt > novo.sam
samtools view -bS novo.sam | samtools sort - reads1
For paired-end or mate-pair reads adjust the parameters accordingly. If we have multiple alignment files these all need to be sorted individually and merged with either samtools or Picard.
Remove PCR duplicates
It is highly recommended to remove PCR duplicates before proceeding with the SNP and Indel calls. We use samtools as an example to remove PCR duplicates.
#remove PCR duplicates from single-end reads
samtools rmdup -s reads1.bam reads1.sorted.bam
#For Paired reads
samtools rmdup reads1.bam reads1.pe.sorted.bam
SNP and Microindel calling
In this example we use the samtools pileup SNP caller and accept reads with a minimum mapping quality of 20:
#The hg18.fasta reference genome sequence is required. samtools.pl is also required
samtools view -u -q 20 reads1.bam | samtools pileup -vc -f hg18.fasta - | samtools.pl varFilter -D 1000 -G 50 > variations.list
Filtering and final list generation
A simple awk script can be used to filter our list of variations where the indel SNP quality is greater than 50 or 60 respectively:
awk '($3=="*"&&$6>=50)||($3!="*"&&$6>30)' variations.list > variations.filtered.list
Exporting to other formats
The samtools distribution contains a helper script called sam2vcf.pl that may be used to convert to the Variant Call Format (VCF):
sam2vcf.pl variations.filtered.list -r hg18.fasta > variations.vcf
Human (hg19, Feb_2009, GRCh37)
single-end, HiSeq2000, Rong Mao Lab
Human (hg18, Mar_2006, NCBI Build 36.1)
DNA-methylation, pair-end, Hiseq2000, Brad Cairns Lab
Chip-seq, single-end, GAIIx, Brad Cairns Lab
C_elegans (CE6, May_2008, WS190)
Small RNA sequencing, single-end, GAIIx, Brenda Bass Lab
A_Thaliana (TAIR8, Mar_2008)
Single-end, HiSeq2000, Jason Stajich Lab
M_musculus (mm9, Jul_2007, NCBI Build 37)
ChIP-Seq, single-end, HiSeq2000, Anne Moon Lab
DNA-methylation, pair-end, HiSeq2000, Brad Cairns Lab
D_rerio (Zv8, Jun_2008)
ChIP-Seq, single-end, HiSeq2000, Brad Cairns Lab
#!/home/ying/py2.6.5/bin/python
"""
Use rsync to synchronize the DATA_HOME(hci-bio2, which raw data was stored) and COMP_HOME(chpc-node,
which process raw data and return output back to source. All the rsync operations are initiated from chpc-node.
Functions:
1. query: query the source if the new data is available
2. download: download new data ( from source to destination)
3. upload: upload result(output) to source.
4. update: update the source by appending a log file from the destination to show the alignment progress on chpc. "begin, percent-finished, end, error, etc"
"""
import os
import time
import operator
import subprocess
DATA_HOME = 'ying@123.123.123.123:~/alignment/'
COMP_HOME = '/home/ying/alignment/'
RSYNC_QUERY_PARAMS = '-are ssh'
RYSNC_DOWNLOAD_PARAMS = '-avuze ssh' #archive, verbose, update-only, compress-during-transfer and trasnfer-over-ssh
RYSNC_UPLOAD_PARAMS = '-avuze ssh' #archive, verbose, update-only, compress-during-transfer and trasnfer-over-ssh
LOG_FILE_NAME = 'log.txt'
def query():
#query_cmd = """rsync -are ssh ying@155.100.235.73:~/alignment/ | awk '{print $1 " " $3"#"$4 " " $5}'"""
cmd = ["rsync",RSYNC_QUERY_PARAMS,DATA_HOME,'|',"""awk '{print $1 " " $3"#"$4 " " $5}'"""]
query_cmd = ' '.join(cmd)
#print query_cmd
p = subprocess.Popen(query_cmd,shell=True,stdout=subprocess.PIPE)
line = p.communicate()
fnames = ''.join(line[:-1]).strip().split('\n')
df = {}
dt = {}
for fname in fnames:
prop,ctime,fn = fname.split(' ')
if fn != '.' and fn != '..':
if prop[0] == 'd': #it is a directory
pass
if prop[0] == '-': #it is a file
parent_path = os.path.dirname(fn)
df.setdefault(parent_path,[]).append(os.path.basename(fn))
#Transform time to seconds, for comparing which job should be run firstly. First In First Run.
#'2010/08/20#16:27:40' -> 1282343260.0
dt[parent_path] = int(time.mktime(time.strptime(ctime,'%Y/%m/%d#%H:%M:%S')))
jobs = []
for k,v in df.items():
if "begin" in v: #the data is ready
if not LOG_FILE_NAME in v: #this job is not running
jobs.append((k,dt[k]))
#print jobs
return sorted(jobs,key=operator.itemgetter(1))
def download(path_to_download):
"""download a whole directory from DATA_HOME to COMP_HOME
"""
normalized_path = path_to_download.strip('/') #/A, /A/ or A/ -> A
source = os.path.join(DATA_HOME,normalized_path)
cmd = ["rsync", RYSNC_DOWNLOAD_PARAMS, source, COMP_HOME] #download the whole directory to COMP_HOME.
download_cmd = ' '.join(cmd)
#print download_cmd
p = subprocess.Popen(download_cmd,shell=True,stdout=subprocess.PIPE)
line = p.communicate()
stdout_message = ''.join(line[:-1]).strip().split('\n')
#print 'Message:',stdout_message
update(normalized_path,stdout_message)
def upload(path,file_to_upload):
"""upload a file from COMP_HOME to DATA_HOME
upload 'file_to_upload' at local path 'folder'.
:/folder/file_to_upload will be uploaded to :/
"""
normalized_path = path.strip('/')+"/" #/A, /A/ or A/ -> A/
source = os.path.join(COMP_HOME,normalized_path,file_to_upload)
dest = os.path.join(DATA_HOME,normalized_path)
cmd = ["rsync",RYSNC_UPLOAD_PARAMS, source,dest]
upload_cmd = ' '.join(cmd)
p = subprocess.Popen(upload_cmd,shell=True,stdout=subprocess.PIPE)
line = p.communicate()
stdout_message = ''.join(line[:-1]).strip().split('\n')
#fn_log = log(os.path.join(COMP_HOME,normalized_src),stdout_message)
update(normalized_path,stdout_message)
def update(path,message):
fn = os.path.join(COMP_HOME,path,LOG_FILE_NAME)
ofn = None
try:
if os.path.exists(fn):
ofn = open(fn,'w+')#append
else:
ofn = open(fn,'w') #new file
print >>ofn,'-'*50
print >>ofn,time.asctime()
for m in message:
print >>ofn,M
ofn.close()
normalized_path = path.strip('/')+"/" #/A, /A/ or A/ -> A/
source = fn
dest = os.path.join(DATA_HOME,normalized_path)
cmd = ["rsync",RYSNC_UPLOAD_PARAMS, source,dest]
subprocess.Popen(cmd,shell=False,stdout=subprocess.PIPE)
except:
pass
def clear():
pass
#print query()
#download('A')
upload('A','hello.txt')
