Charpter 15 Advanced Command Line
背景
在Linux下的數(shù)據(jù)分析從質(zhì)控到比對到后續(xù)分析等猩谊，每一步都需要獨立的軟件去實現(xiàn)，并且在一些步驟之間還需要對數(shù)據(jù)的格式進行轉(zhuǎn)換。而想要整合這些軟件去實現(xiàn)“一鍵式”的分析，就需要去學習一些進階腳本的寫法了。
編程語言包括兩類：編譯型語言(complied programs)和解釋型語言(interpreted programs)。編譯型語言更符合底層計算機語言裙盾，故程序運行效率高速度快，但上手存在一定的困難他嫡。解釋型語言更符合人類語言番官，有其特殊的解釋器來解析程序語言，上手較方便钢属，但程序運行效率沒有編譯型語言高徘熔。
生物信息里生物汪首先學習的還是編譯型語言：主要包括Perl，Python淆党，R酷师，AWK等∪疚冢基本語言運用熟練了有需求再看看編譯型語言C山孔，能更好的理解計算機數(shù)據(jù)處理。
AWK
入門時學習Perl的慕匠，也練習下用perl單行實現(xiàn)書中AWK的效果吧
### 數(shù)據(jù)下載
wget http://data.biostarhandbook.com/bam/demo.bam
samtools view demo.bam |cut -f 1,3,4 |head
### AWK看指定行饱须，并進行簡單的計算
samtools view demo.bam | awk '{ print $3,$2,$1, 10^(-$5/10) }' |head
###perl 單行實現(xiàn)
samtools view demo.bam | perl -alne ' print "$F[2]\t$F[1]\t$F[0]",10**(-$F[4]/10) '
AWK 基本格式 :awk 'CONDITION { ACTIONS }', awk 'CONDITION1 { ACTIONS1 } CONDITION2 { ACTIONS2 } CONDITION3 { ACTIONS3 }'
### AWK 篩選指定列小于60的值
samtools view demo.bam | awk '$5 < 60 { print $5 }'
##Perl oneliner
samtools view demo.bam | perl -alne 'print $F[4] if $F[4] < 60'
AWK默認是對每行的空白部分(spaces, tabs)為分隔符
- 指定分隔符
awk -F '\t'，相同與Perl單行的-F(需與-a參數(shù)一起)台谊。 - NR 指行數(shù)，number of row
- NF 指列數(shù)譬挚， number of columns
###計算TLEN插入片段長度大于0的平均值锅铅。
samtools view -f 2 demo.bam |awk ' $9 > 0 {sum+=$9 ; count+=1} END{print sum/count} '
### perl oneliner
samtools view -f 2 demo.bam |perl -alne ' if($F[8]>0){$sum+=$F[8];$count+=1 } END{print $sum/$count}'
## 加了個time 看awk比perl運行速度更快一些
BioAwk：-c參數(shù)等
##檢測sam中cigar有deletion
samtools view -f 2 demo.bam |bioawk -c sam '$cigar ~ /D/{print $cigar}'|wc -l
## perl 單行
samtools view -f 2 demo.bam | perl -alne ' print $F[5] if $F[5]=~/D/ ' |wc -l
不同condition和不同的action：
samtools depth demo.bam | awk ' $3 <= 3 { $1="LO" } $3 > 3 { $1="HI" } { print $0 }'
## perl oneliner
samtools depth demo.bam |perl -alne 'if($F[2]<=3){$F[0]="LO";print "$F[0]\t$F[1]\t$F[2]"}else{$F[0]="HI"; print "$F[0]\t$F[1]\t$F[2]"}'
腳本scripts
bash腳本的一些基本用法：
- 開頭
#!/bin/bash - 設置糾錯參數(shù)(Dead programs tell no lies):
set -ueo pipefail - 基本變量
NAME=JOHN.... ${NAME}(注意變量賦值=前后不能有空格,使用變量用${NAME}) - 命令行參數(shù)：
NAME=$1(相當于perl里的$ARGV[0]) -
對于絕對地址的文件：
- 文件名(basename{FULL})
- 后綴${FULL#*.}
- 前綴${FULL%.*}
- 將命令結(jié)果賦予一個變量`ls -1 |wc -l`
-
ls -1-1參數(shù)的使用 -
并行處理parallel：
parallel --verbose --eta- --verbose : 指具體并行處理的命令行
- --eta：指具體處理的時間
find . -name '*.fastq' -print0 |parallel -0 --verbose "grep ATGC {} >{}.out"ls * |parallel -i -verbose "grep ATGC {} > {}.OUT"
### bash的一些詭異變量處理
# Suppose this is the full path.
FULL=/data/foo/genome.fasta.tar.gz
# To make it: genome.fasta.tar.gz
NAME=$(basename ${FULL})
# To make it: fasta.tar.gz
EXT1=${FULL#*.}
# To get only the extension: gz
EXT2=${FULL##*.}
# To get the other side: /data/foo/genome.fasta.tar
START=${FULL%.*}
START=${FULL%%.*}
### 獲取SRR數(shù)據(jù)10000行并質(zhì)控
#
# Usage: getproject.sh PRJN NUM
#
# Example: bash getproject.sh PRJNA257197 3
#
# This program downloads a NUM number of experiments from an SRA Bioproject
# identified by the PRJN number then runs FastQC quality reports before
# and after trimming the data for quality.
# Immediately stop on errors.
set -ueo pipefail
# The required parameter is the run id.
PRJN=$1
# How many experimental runs to get.
NUM=$2
# What is the conversion limit for each data.
LIMIT=10000
# This is an internal variable that holds the selected ids.
SRA=short.ids
# Keep only the required number of ids.
cat ids.csv | head -${NUM} > $SRA
# Place an echo before each command to see what will get executed when it runs.
cat $SRA | xargs -n 1 -I {} echo fastq-dump --split-files -X ${LIMIT} {}
# Generate the commands for fastqc.
cat $SRA | xargs -n 1 -I {} echo fastqc {}_1.fastq {}_2.fastq
# Generate the commands for trimmomatic.
# Here we run it with the -baseout flag to generatevfile extension of .fq.
cat $SRA | xargs -n 1 -I {} echo trimmomatic PE -baseout {}.fq {}_1.fastq {}_2.fastq SLIDINGWINDOW:4:30
# Run FastQC on the new data that matches the *.fq extension.
cat $SRA | xargs -n 1 -I {} echo fastqc {}_1P.fq {}_2P.fq
Charpter 16 DATA Visualization
背景
永遠不能輕易的就相信電腦，大腦才是最強有力的分析工具减宣。
能夠可視化的數(shù)據(jù)包括：FASTA, BED, GFF, SAM/BAM
常用的數(shù)據(jù)可視化工具包括：
- IGV(Intergrative Genomics Viewer)
- SeqMonk 用來分析RNA-seq比對數(shù)據(jù)的工具
- UCSC Genome Browser
- Ensembl Genome Browser
