Biostar_handbook||charpter 15_16

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