Linux题库解答

前言

这个教程源于曾健明，题库地址为其博客。所使用的Linux系统是曾健明的服务器，服务器信息如下：

/usr/local/bin/miniconda3/bin路径下面安装了生物信息学软件，可以使用全路径调用它们，或者添加该目录到环境变量。
两个练手数据：
2.1 转录组数据：/public/study/mRNAseq/tair/的转录组的测试数据，具体教程可以看其博客：http://www.bio-info-trainee.com/2809.html，可以用ls -lh --color=auto /public/study/mRNAseq/tair/ 命令查看，然后所有的处理，在自己的目录处理，不需要拷贝那些数据.
2.2 找变异：/public/study/variant-calling/tair。
可以参见教程：http://www.biotrainee.com/thread-696-1-1.html

其他的NGS流程，参见Github：https://github.com/jmzeng1314/NGS-pipeline。

基础知识

文件目录操作

系统管理

1.查看当前登录用户w
mark

2.查看命令历史l~/.bash_history
Bash shell在“~/.bash_history”（“~/”表示用户目录）文件中保存了500条使用过的命令，这样可以使你输入使用过的长命令变得容易。每个在系统中拥有账号的用户在他的目录下都有一个“.bash_history”文件。

用户权限

文本操作

题库测试

1. 在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列。

mkdir -p 1/2/3/4/5/6/7/8/9
mkdir是建立文件夹命名，添加-p参数可以直接创建多级文件夹，如下所示：
mark

2.在创建好的文件夹下面，比如我的是 /Users/jimmy/tmp/1/2/3/4/5/6/7/8/9 ，里面创建文本文件 me.txt

先切换到/1/2/3/4/5/6/7/8/9，直接输入cd，然后不停地按Tab键，就会进入下一层文件夹。
touch me.txt创建txt文件。
mark

3.在文本文件 me.txt 里面输入内容:
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?

输入vim me.txt，输入上述文字，然后按Esc键，输入wp退出。

4. 删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt

rm -rf *

5. 在任意文件夹下面创建 folder1~5这5个文件夹，然后每个文件夹下面继续创建 folder1~5这5个文件夹，效果如下：

mkdir -p folder_{1..5}/folder_{1..5}
如下所示：
mark

6. 在第五题创建的每一个文件夹下面都创建第二题文本文件 me.txt ，内容也要一样。

先创建一个tmp3目录，其目录下有5个文件夹，5个文件夹下面又有5个文件夹，然后将me.txt拷贝到这25个文件夹中：

1 2	mkdir -p tmp3/folder_{1..5}/folder_{1..5} echo tmp3/folder_{1..5}/folder_{1..5} \| xargs -n 1 cp -v me.txt

mark

7. 再次删除掉前面几个步骤建立的文件夹及文件。
rm -rf *

8. 下载 http://www.biotrainee.com/jmzeng/igv/test.bed 文件，后在里面选择含有 H3K4me3 的那一行是第几行，该文件总共有几行。

1
2
3

wget http://www.biotrainee.com/jmzeng/igv/test.bed
grep -n "H3K4me3" test.bed
wc test.bed

结果如下：
mark
可知，H3K4me3在第8行，用Notepad++打开的话，核对一下，确实如此。
mark
用wc命令查看有多少行，一共10行。

9. 下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件，并且解压，查看里面的文件夹结构

1
2
3

wget http://www.biotrainee.com/jmzeng/rmDuplicate.zip
unzip rmDuplicate.zip
cd rmDuplicate.zip

查看文件：
mark

10. 打开第九题解压的文件，进入 rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚生物信息学里面的SAM/BAM 定义是什么。
mark

SAM文件与BAM文件

SAM是一种序列比对格式标准，由sanger制定，是以TAB为分割符的文本格式。主要应用于测序序列mapping到基因组上的结果表示，当然也可以表示任意的多重比对结果。

SAM要处理好的问题：

非常多序列（read)，mapping到多个参考基因组（reference）上；
同一条序列，分多段（segment）比对到参考基因组上；
无限量的，结构化信息表示，包括错配、删除、插入等比对信息；
SAM分为两部分，注释信息（header section）和比对结果部分（alignment section），注释信息可有可无，都是以@开头，用不同的tag表示不同的信息，主要有@HD，说明符合标准的版本、对比序列的排列顺序；@SQ，参考序列说明；@RG，比对上的序列（read）说明；@PG，使用的程序说明；@CO，任意的说明信息。

mark

>
比对结果部分（alignment section），每一行表示一个片段（segment）的比对信息，包括11个必须的字段（mandatory fields）和一个可选的字段，字段之间用tag分割。必须的字段有11个，顺序固定，不可用时，根据字段定义，可以为’0‘或者’*‘，这是11个字段包括：

QNAME，比对片段的（template）的编号；
FLAG，位标识，template mapping情况的数字表示，每一个数字代表一种比对情况，这里的值是符合情况的数字相加总和；
RNAME，参考序列的编号，如果注释中对SQ-SN进行了定义，这里必须和其保持一致，另外对于没有mapping上的序列，这里是’‘；
POS，比对上的位置，注意是从1开始计数，没有比对上，此处为0；
MAPQ，mappint的质量；
CIGAR，简要比对信息表达式（Compact Idiosyncratic Gapped Alignment Report），其以参考序列为基础，使用数字加字母表示比对结果，比如3S6M1P1I4M，前三个碱基被剪切去除了，然后6个比对上了，然后打开了一个缺口，有一个碱基插入，最后是4个比对上了，是按照顺序的；
RNEXT，下一个片段比对上的参考序列的编号，没有另外的片段，这里是’‘，同一个片段，用’=‘；
PNEXT，下一个片段比对上的位置，如果不可用，此处为0；
TLEN，Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0；
SEQ，序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；
QUAL，序列的质量信息，格式同FASTQ一样。
可选字段（optional fields)，格式如：TAG:TYPE:VALUE，其中TAG有两个大写字母组成，每个TAG代表一类信息，每一行一个TAG只能出现一次，TYPE表示TAG对应值的类型，可以是字符串、整数、字节、数组等。

11. 安装 samtools 软件

mkdir soft
cd soft
wget https://github.com/samtools/samtools/releases/download/1.7/samtools-1.7.tar.bz2
tar jxvf samtools-1.7.tar.bz2
cd samtools-1.7/
./configure --prefix=/trainee/home/jlzhang/biosoft
make
make install
export PATH=$PATH:/trainee/home/jlzhang/biosoft/bin:$PATH >>~/.bashrc
source ~/.bashrc

可以通过echo $PATH来查看用户变量。

12. 打开后缀为BAM 的文件，找到产生该文件的命令。

用samtools view打开BAM文件，过程如下：

进入相应文件夹，即samtools view tmp.rmdup.bam | head
输入命令：samtools view tmp.rmdup.bam | head
打开后的BAM文件如下所示：
mark

13. 根据上面的命令，找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体。
从网上找到的信息，hg38含有24条染色体，统计代码如下所示：

先找到hg38.fa文件，如下所示：

1	locate hg38.fa

mark
随便挑一个进行统计，如下所示：

1 2	grep "chr[1-9,X,Y]" /umac/reference/genome/hg38/hg38.fa\| cut -d_ -f1 > statistic.chrom sort statistic.chrom \| uniq \| wc

可以发现，hg38有24条染色体，如下所示：
mark

其中,sort是对文本文件中的字符进行排序，uniq是去掉重复的行，但是uniq只能去掉相邻的重复行，因此之前要进行sort。

14. 上面的后缀为BAM 的文件的第二列，只有 0 和 16 两个数字，用 cut/sort/uniq等命令统计它们的个数。

samtools view tmp.rmdup.bam > bamfile
awk '{print $2}' bamfile > zero_one
grep -o "0" zero_one |wc -l
grep -o "16" zero_one |wc -l

mark

15. 重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件，再次查看第二列，并且统计

1
2
3

samtools view tmp.rmdup.bam > pairbam
awk '{print $2}' pairbam > number 
awk '{print $2}' pairbam |sort|uniq -c

mark

16. 下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，查看里面的文件夹结构，这个文件有2.3M，注意留心下载时间及下载速度。

wget http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
ls -lh
unzip sickle-results.zip
cd sickle-results
ls -lh

mark

17. 解压 sickle-results/single_tmp_fastqc.zip 文件，并且进入解压后的文件夹，找到 fastqc_data.txt 文件，并且搜索该文本文件以 >>开头的有多少行？

cd sickle-results
unzip single_tmp_fastqc.zip
cd single_tmp_fastqc
grep ">>"  fastqc_data.txt |wc -l

一共是24行
mark

18. 下载 http://www.biotrainee.com/jmzeng/tmp/hg38.tss 文件，去NCBI找到TP53/BRCA1等自己感兴趣的基因对应的 refseq数据库 ID，然后找到它们的hg38.tss 文件的哪一行。
https://www.ncbi.nlm.nih.gov/gene/7157

打开PUBMED，选GENE，输入“TP53”，点击Search，如下所示：

Pubmed中有很多TP53的信息，但是hg38是人类基因组，因此选择第一行的TP53就行，进入即可，找到Go to reference sequence details，点进去即可，如下所示：
mark

找到Go to reference sequence details，进去，其中有NM开头的编号，NM表示mRNA的RefSeq，TP53的就是NM_00054635（小数点部分是版本号，不用管），在hg38.tss中查询即可，BRCA1按照同样操作，查询结果如下所示：

1
2
3

wget http://www.biotrainee.com/jmzeng/tmp/hg38.tss
grep "NM_000546" hg38.tss 
grep "NM_005905" hg38.tss

mark

19. 解析hg38.tss 文件，统计每条染色体的基因个数。

hg38.tss文件是基因的坐标，查看了一下文件，文件内容如下所示：
mark

第1列是基因名称，第2列是染色体名称，第3列到第4列是转录起止坐标，第5列不清楚。

1
2
3

awk '$5 > 0' hg38.tss > gene2
awk '{print $2}' gene2 > gene
cut -d_ -f1 gene |sort |uniq -c

mark

20. 解析hg38.tss 文件，统计NM和NR开头的序列，了解NM和NR开头的含义。

没在hg38.tss文件里找到相应的序列，只找到了名称，提取NM和NR开头的名称命令如下：

1 2	grep "NM" hg38.tss grep "NR" hg38.tss

NM，mRNA mixed，转录组产物序列；成熟mRNA转录本序列。
NR，RNA mixed，非编码的转录子序列，包括结构RNAs，假基因转子等。

参考资料：