Charpter 15 Advanced Command Line
背景
在Linux下的数据分析从质控到比对到后续分析等，每一步都需要独立的软件去实现，并且在一些步骤之间还需要对数据的格式进行转换。而想要整合这些软件去实现“一键式”的分析，就需要去学习一些进阶脚本的写法了。
编程语言包括两类：编译型语言(complied programs)和解释型语言(interpreted programs)。编译型语言更符合底层计算机语言，故程序运行效率高速度快，但上手存在一定的困难。解释型语言更符合人类语言，有其特殊的解释器来解析程序语言，上手较方便，但程序运行效率没有编译型语言高。
生物信息里生物汪首先学习的还是编译型语言：主要包括Perl，Python，R，AWK等。基本语言运用熟练了有需求再看看编译型语言C，能更好的理解计算机数据处理。
AWK
入门时学习Perl的，也练习下用perl单行实现书中AWK的效果吧
### 数据下载
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 单行实现
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参数一起)。 - NR 指行数，number of row
- NF 指列数， 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参数等
##检测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 - 设置纠错参数(Dead programs tell no lies):
set -ueo pipefail - 基本变量
NAME=JOHN.... ${NAME}(注意变量赋值=前后不能有空格,使用变量用${NAME}) - 命令行参数：
NAME=$1(相当于perl里的$ARGV[0]) -
对于绝对地址的文件：
- 文件名(basename{FULL})
- 后缀${FULL#*.}
- 前缀${FULL%.*}
- 将命令结果赋予一个变量`ls -1 |wc -l`
-
ls -1-1参数的使用 -
并行处理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数据10000行并质控
#
# 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
背景
永远不能轻易的就相信电脑，大脑才是最强有力的分析工具。
能够可视化的数据包括：FASTA, BED, GFF, SAM/BAM
常用的数据可视化工具包括：
- IGV(Intergrative Genomics Viewer)
- SeqMonk 用来分析RNA-seq比对数据的工具
- UCSC Genome Browser
- Ensembl Genome Browser
