옥탑방람보

< 선행설치>

1. autoconf 설치 (2.65 이상의 버전)

다운로드 - http://ftp.gnu.org/gnu/autoconf/

설치 - ./configure, make, make install

(필요할경우) vi /etc/profile 에서 export PATH=/usr/local/bin:$PATH 추가

2. automake 설치

다운로드 - http://ftp.gnu.org/gnu/automake/

설치 - ./configure, make, make install

3. m4 설치

다운로드 - http://ftp.gnu.org/gnu/m4/

설치 - ./configure, make, make install

4. SSE2 지원 확인

cat /proc/cpuinfo 시에 flags 부분에 sse2 가 있는지 확인

5. git 설치 및 TMAP 다운로드

$> yum install git-XXX

$> git clone git://github.com/iontorrent/TMAP.git  (방화벽 9418 오픈되어 있어야함)

$> cd TMAP

$> git submodule init

$> git submodule update

<TMAP 설치>

1. 설치

$> cd TMAP

$> chmod +x autogen.sh

$> chown -R user.group ../TMAP

$root> sudo ./autogen.sh (perl모듈 사용으로 root로 해야하는 경우 있음)

$> mkdir ~/install/TMAP_130122

$>./configure --prefix=~/install/TMAP_130122

$root> sudo make (perl모듈 사용으로 root로 해야하는 경우 있음)

$> make install

reference site: http://samtools.sourceforge.net/
manual : http://samtools.sourceforge.net/SAM1.pdf

The Samtools is an essential tool to handle files with sam (bam) format.

99 (decimal) -> 1100011 (binary) : paired, proper pair, mapped, forward, mate reverse .... (Watch the number one by one by reverse order)

FROM: http://sourceforge.net/apps/mediawiki/samtools/index.php?title=SAM_FAQ

About SAMtools

How to convert SAM to BAM?
If your SAM file has header @SQ lines, you may get BAM by:
samtools view -bS aln.sam > aln.bam
If not, you need to have your reference file ref.fa and then do this:
samtools faidx ref.fa
samtools view -bt ref.fa.fai aln.sam > aln.bam
The second method also works if your SAM file has @SQ lines. After conversion, you would probably like to sort and index the alignment to enable fast random access:
samtools sort aln.bam aln-sorted
samtools index aln-sorted.bam

I want to call SNPs and short indels.
For a short answer, do this:
samtools pileup -vcf ref.fa aln.bam | tee raw.txt | samtools.pl varFilter -D100 > flt.txt
awk '($3=="*" && $6>=50)||($3!="*" && $6>=20)' flt.txt > final.txt
For a long answer, see SAM protocol. Please always remember to set the maximum depth (-D) in filtering.

I want to call SNPs from one chromosome only.
Index your alignment with the `index' command and:
samtools view -u aln.bam chr10 | samtools pileup -vcf ref.fa - > chr10.raw.txt
Please read [http://samtools.sourceforge.net/pipe.shtml this page] for more information.

The integer FLAG field is not friendly to eyes.
You may get string FLAG by:
samtools view -X aln.bam | less -S
For more information, please check out:
samtools view -?


Pileup output.

This is explained in the [http://samtools.sourceforge.net/samtools.shtml manual page]. Or briefly (when you invoke pileup with the -c option):
1. reference sequence name
2. reference coordinate
3. reference base, or `*' for an indel line
4. genotype where heterozygotes are encoded in the [http://biocorp.ca/IUB.php IUB code]: M=A/C, R=A/G, W=A/T, S=C/G, Y=C/T and K=G/T; indels are indicated by, for example, */+A, -A/* or +CC/-C. There is no difference between */+A or +A/*.
5. Phred-scaled likelihood that the genotype is wrong, which is also called `consensus quality'.
6. Phred-scaled likelihood that the genotype is identical to the reference, which is also called `SNP quality'. Suppose the reference base is A and in alignment we see 17 G and 3 A. We will get a low consensus quality because it is difficult to distinguish an A/G heterozygote from a G/G homozygote. We will get a high SNP quality, though, because the evidence of a SNP is very strong.
7. [http://en.wikipedia.org/wiki/Root_mean_square root mean square] (RMS) mapping quality
8. # reads covering the position
9. read bases at a SNP line (check the manual page for more information); the 1st indel allele otherwise
10. base quality at a SNP line; the 2nd indel allele otherwise
11. indel line only: # reads directly supporting the 1st indel allele
12. indel line only: # reads directly supporting the 2nd indel allele
13. indel line only: # reads supporting a third indel allele
If pileup is invoked without `-c', indel lines and columns between 3 and 7 inclusive will not be outputted.

I see `*' in the pileup sequence column. What are they?
A star at the sequence column represents a deletion. It is a place holder to make sure the number of bases at that column matches the read depth column. Simply ignore `*' if you do not use this information.


Does samtools generate the consensus sequence like Maq?

Yes. Try this:
samtools pileup -cf ref.fa aln.bam | samtools.pl pileup2fq -D100 > cns.fastq
Again, remember to set -D according to your read depth. Note that pileup2fq applies fewer filters in comparison to varFilter, and you may see tiny inconsistency between the two outputs.

I want to get `unique' alignments from SAM/BAM.
We prefer to say an alignment is `reliable' rather than `unique' as `uniqueness' is not well defined in general cases. You can get reliable alignments by setting a threshold on mapping quality:
samtools view -bq 1 aln.bam > aln-reliable.bam
You may want to set a more stringent threshold to get more reliable alignments.

In repetitive regions, SAMtools call all bases as 'A' although there are no 'A' bases in reads.
This is due to a floating underflow in the MAQ SNP calling model used by default and only happens in repetitive regions. These calls are always filtered out. However, if you are uncomfortable with this, you may use the simplified SOAPsnp model with:
samtools -avcf ref.fa aln.bam > raw.txt
The MAQ model and SOAPsnp model usually deliver very similar SNP calls.

How are SNPs and indels called and filtered by SAMtools?
By default, SNPs are called with a Bayesian model identical to the one used in MAQ. A simplified SOAPsnp model is implemented, too. Indels are called with a simple Bayesian model. The caller does local realignment to recover indels that occur at the end of a read but appear to be contiguous mismatches. For an example, see [http://samtools.sourceforge.net/images/seq2-156.png this picture].

The varFilter filters SNPs/indels in the following order:
* d: low depth
* D: high depth
* W: too many SNPs in a window (SNP only)
* G: close to a high-quality indel (SNP only)
* Q: low root-mean-square (RMS) mapping quality (SNP only)
* g: close to another indel with more evidence (indel only)
The first letter indicates the reason why SNPs/indels are filtered when you invoke varFilter with the `-p' option. A SNP/indel filtered by a rule higher in the list will not be tested against other rules.

• AA    ancestral allele
• AC    allele count in genotypes, for each ALT allele, in the same order as listed
• AF    allele frequency for each ALT allele in the same order as listed: use this when estimated from primary data, not called genotypes
• AN    total number of alleles in called genotypes
• BQ    RMS base quality at this position
• CIGAR    cigar string describing how to align an alternate allele to the reference allele
• DB    dbSNP membership
• DP    combined depth across samples, e.g. DP=154
• END    end position of the variant described in this record (esp. for CNVs)
• H2    membership in hapmap2
• MQ    RMS mapping quality, e.g. MQ=52
• MQ0    Number of MAPQ == 0 reads covering this record
• NS    Number of samples with data
• SB    strand bias at this position
• SOMATIC    indicates that the record is a somatic mutation, for cancer genomics
• VALIDATED    validated by follow-up experiment

• I16
  • 1 #reference Q13 bases on the forward strand 2 #reference Q13 bases on the reverse strand
  • 3 #non-ref Q13 bases on the forward strand 4 #non-ref Q13 bases on the reverse strand
  • 5 sum of reference base qualities 6 sum of squares of reference base qualities
  • 7 sum of non-ref base qualities 8 sum of squares of non-ref base qualities
  • 9 sum of ref mapping qualities 10 sum of squares of ref mapping qualities
  • 11 sum of non-ref mapping qualities 12 sum of squares of non-ref mapping qualities
  • 13 sum of tail distance for ref bases 14 sum of squares of tail distance for ref bases
  • 15 sum of tail distance for non-ref bases 16 sum of squares of tail distance for non-ref
• INDEL    Indicating the variant is an INDEL.
• DP    The number of reads covering or bridging POS.
• DP4    Number of 1) forward ref alleles; 2) reverse ref; 3) forward non-ref; 4) reverse non-ref alleles, used in variant calling. Sum can be smaller than DP because low-quality bases are not counted.
• PV4    P-values for 1) strand bias (exact test); 2) baseQ bias (t-test); 3) mapQ bias (t); 4) tail distance bias (t)
• FQ    Consensus quality. If positive, FQ equals the phred-scaled probability of there being two or more different alleles. If negative, FQ equals the minus phred-scaled probability of all chromosomes being identical. Notably, given one sample, FQ is positive at hets and negative at homs.
• AF1    EM estimate of the site allele frequency of the strongest non-reference allele.
• CI95    Equal-tail (Bayesian) credible interval of the site allele frequency at the 95% level.
• PC2    Phred-scaled probability of the alternate allele frequency of group1 samples being larger (,smaller) than of group2 samples.
• PCHI2    Posterior weighted chi^2 P-value between group1 and group2 samples. This P-value is conservative.
• QCHI2    Phred-scaled PCHI2
• RP    Number of permutations yeilding a smaller PCHI2