Linux 生信项目案例在线分析
本文转自https://freeaihub.com/article/bioinformatics-in-linux.html，该页可在线进行操作。
Linux对生物信息的学习和实践有强大的辅助作用，不管自己编写shell命令脚本处理数据，还是使用现有丰富生物信息分析工具。Linux强大命令行功能，可以快速、批量、灵活的处理数据的提取、统计和整理等耗时耗力的重复性工作。日常生信分析中，多数整理工作都是用Linux命令的组合完成的，这相比于写完整的Python或Perl程序更简便快捷。
本节我们将通过一个案例来了解如何使用Linux来进行Linux生项项目的实践分析。
bowtie2简介与在线安装
Bowtie2 是将测序reads与长参考序列比对工具 (适用于将长度大约为50到100或1000字符的reads与相对较长的基因组, 如哺乳动物,进行比对)。
通常是比较基因组学(包括识别变体（variation calling），ChIP-seq，RNA-seq，BS-seq)管道的第一步。
可以处理非常长的读数（即10s或100s的千字节），但它针对近期测序仪产生的读数长度和误差模式进行了优化，如Illumina HiSeq 2000，Roche 454和Ion Torrent仪器。
Bowtie2使用FM索引（基于Burrows-Wheeler Transform 或 BWT）对基因组进行索引，以此来保持其占用较小内存。
#安装bowtie2
cd /usr/local
cp /share/tar/bowtie2-2.3.4.1-linux-x86_64.zip .
unzip bowtie2-2.3.4.1-linux-x86_64.zip
mv bowtie2-2.3.4.1-linux-x86_64 bowtie2
echo 'export PATH=/usr/local/bowtie2:$PATH' >> ~/.bashrc
source ~/.bashrc #生效环境变量
#验证安装是否成功
bowtie2 --version
samtool简介与在线安装
samtools是一个用于操作sam和bam文件的工具合集。能够实现二进制查看、格式转换、排序及合并等功能，结合sam格式中的flag、tag等信息，还可以完成比对结果的统计汇总。同时利用linux中的grep、awk等操作命令，还可以大大扩展samtools的使用范围与功能。包含有许多命令。
#安装samtool
cd /usr/local
cp /share/tar/samtools-1.9.tar.bz2 .
#wget -c https://github.com/samtools/samtools/releases/download/1.9/samtools-1.9.tar.bz2
tar xvf samtools-1.9.tar.bz2
/usr/local/samtools-1.9/configure --prefix=/usr/local/samtools-1.9
echo 'export PATH=/usr/local/samtools-1.9:$PATH' >> ~/.bashrc
source ~/.bashrc #生效环境变量
samtools --version
下载相应分析数据
下载http://www.biotrainee.com/jmzeng/igv/test.bed 文件，后在里面选择含有 H3K4me3 的那一行是第几行，该文件总共有几行。
pwd
wget –c http://www.biotrainee.com/jmzeng/igv/test.bed
ls
#(-n 标记行数，-o 只显示匹配上的，--color匹配文字出现颜色）
grep -n -o --color H3K4me3 test.bed
#(wc显示文件的行数、单词数、字节数）
cat test.bed |wc -l
下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件，并且解压，查看里面的文件夹结构
wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip
ls
unzip rmDuplicate.zip
ls
tree rmDuplicate
打开解压的文件，进入 rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚生物信息学里面的SAM/BAM定义是什么【另外关于sam头部注释信息在tmp.header中】
SAM的全称是sequence alignment/map format, 主要用于测序序列mapping到基因组上的结果表示
BAM是SAM的二进制压缩文件，不能像sam可以用less查看，它需要用samtools view打开
SAM分数据通常有两部分，头部注释信息（header section ）和主体比对结果部分 （alignment section）
NM:i 经过编辑的序列长度(edit distance)
MD:Z 错配位置/碱基(mismatching positions/bases)
AS:i 匹配得分(Alignment score)
XS:i 第二好的匹配得分(suboptimal alignment score)
YS:i mate序列匹配的得分
BC: 条码序列（barcode sequence)
rmDuplicate/samtools/single中的sorted.bam文件中，统计第二列各flag个数
cd ~/rmDuplicate/samtools/single
samtools view tmp.sorted.bam | cut -f2 | sort -n | uniq -c
下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，进入解压后的文件夹test1_fastqc，找到 fastqc_data.txt 文件，统计 以>>开头的有多少行
cd ~
wget http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
unzip sickle-results.zip
cd sickle-results
unzip test1_fastqc.zip
cd test1_fastqc
cat fastqc_data.txt | grep '>>' | wc -l
下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss (存储转录起始位点信息)文件，在NCBI中找到自己感兴趣的基因对应的 refseq数据库 ID（这里我以乳腺癌BRCA1为例），然后找到它在hg38.tss 文件的哪一行
cd ~
wget http://www.biotrainee.com/jmzeng/tmp/hg38.tss
cat hg38.tss | grep 'NM_007294'
关于RefSeq数据库：
RefSeq数据库中所有的数据是一个非冗余的、提供参考标准的数据，包括染色体、基因组（细胞器、病毒、质粒）、蛋白、RNA等。refseq序列是NCBI筛选过的非冗余数据库，比GeneBank准确。
关于它的ID：NM开头的表示标准序列，XM表示预测的蛋白编码序列，NR表示非编码蛋白的mRNA序列，AF开头的表示克隆序列，BC开头的表示模板序列
另外，你可能见过gi|4557284|ref|NM_000646.1|[4557284]这种格式
gi就是代表genebank identifier；ref就是对应的refseq中的ID啦
解析hg38.tss 文件，统计每条染色体的基因个数
cd ~
cut -f2 hg38.tss |cut -d'_' -f1 | sort |uniq -c | sort -rn
统计hg38.tss 文件中NM和NR开头的序列，了解NM和NR开头的含义
#统计NM开头
cat hg38.tss | grep 'NM' | wc -l
#同理可以统计NR开头
#NM开头的表示标准序列，可以转录成蛋白质的基因
#NR非编码蛋白的mRNA序列
本节总结了两款生信专用软件bowtie2和samtool的安装，并使用了Linux中的命令，如wc,cut等常见命令完成生物信息项目的简单分析。
参考:1.刘小泽 《Linux学习》 简书 https://www.jianshu.com/p/0b0e7c0f31db
