生信linux20题作业学习笔记

转载自猪猪生信（公众号）生信人应该这样学linux-linux20题作业
https://mp.weixin.qq.com/s/3b5NORO3wWZrE1HV4GQMNA
1.在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列。

mkdir -p 1/2/3/4/5/6/7/8/9 #递归创建文件夹
tree #查看目录结构

2.创建好的文件夹下面，里面创建文本文件 me.txt

cd /1/2/3/4/5/6/7/8/9 #进入目录
vi me.txt 

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

vi me.txt
i #输入i插入后直接输内容
wq #保存退出文本

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

cd #返回自己文件夹所在的初始目录
rm -rf 1/ #批量递归删除目录及文件。

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

mkdir -p test{1..5}/test{1..5}

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

vi test.sh #写一个shell脚本运行
#脚本中内容
 for i in {1..5};do
     cd /root/linux_test_20/folder$i
     for I in {1..5};do
     cd /root/linux_test_20/folder$i/folder$I
     echo -e "Go to: http://www.biotrainee.com/ \n I love bioinfomatics. \n And you ?" > me.txt
     done
  done
sh test.sh #运行脚本

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

rm -rf * #注意所在目录，慎重使用啊

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

wget -c http://www.biotrainee.com/jmzeng/igv/test.bed #可以断点续传下载
grep -no H3K4me3 test.bed #n显示行数 o只显示匹配到的字符
wc -l test.bed #统计行数

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

wget -c  http://www.biotrainee.com/jmzeng/rmDuplicate.zip
unzip rmDuplicate.zip #解压
tree #下载

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

#sam/bam:sam文件是测序的短序列比对后带有比对信息的文件，bam是为了减小文件的大小方便存储，由sam文件转换为二进制的二进制文件
cd rmDuplicate/samtools/single
less -s *.sam

11.安装 samtools 软件

#使用conda安装
#先创建samtools的小环境，再独立的小环境安装是为了避免与其他软件产生冲突
conda create -n samtools
#激活小环境
conda activate samtools
#安装
conda install samtools

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

samtools view -h tmp.rmdup.bam| less

NGS分析中大多数文件都是由header和record两部分组成，加上-h参数后可以将header显示出来，默认是不显示的
header内容
@HD：是必须的标准文件头，包含版本信息；
@SQ：参考序列染色体名字和长度信息 (SN:scaffold name; LN: length)；
@RG：重要read group信息，通常包含测序平台，测序文库和样本ID等信息，分析时用于区分不同样本（重测序时用到）；
@PG：生成此文件的操作过程和参数信息 （program）。

record内容
每一行就是一条read比对上参考基因组的信息，总共12列，用tab键分割。

1. read名称；
2. 比对信息位flag值；
3. 参考序列染色体编号；
4. 5′端起始位置；
5. MAPQ：mapping quality，描述比对的质量，数字越大，特异性越高；
6. CIGAR字符串，记录插入、删除、错配等信息；
7. 配对read所比对到的染色体，仅双端测序的数据才有；
8. 配对read所比对到的位置，仅双端测序的数据才有；
9. 插入片段的长度，仅双端测序的数据才有；
10. read序列；
11. read质量值；
12. 12列以后的信息都是metadata，程序用标记

提示一下命令是：（在PG行）

/home/jianmingzeng/biosoft/bowtie/bowtie2-2.2.9/bowtie2-align-s --wrapper basic-0 -p 20 -x 
/home/jianmingzeng/reference/index/bowtie/hg38 -S 
/home/jianmingzeng/data/public/allMouse/alignment/WT_rep2_Input.sam -U /tmp/41440.unp

十三题、根据上面的命令，找到 /rmDuplicate/samtools/single/tmp.sorted.bam 文件中具体有多少条染色体

samtools view  tmp.sorted.bam | awk {'print $2'} | cut -c4-9 | sort | uniq -c | grep -v "_"

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

samtools view  tmp.sorted.bam | cut -f 2 |sort |uniq -c

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

samtools view samtools/paired/tmp.sorted.bam | cut -f 2 | sort | uniq -c

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

wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip

du -sh sickle-results.zip

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

unzip sickle-results.zip
unzip single_tmp_fastqc.zip
cat fastqc_data.txt | awk '/^>>/{print $0}' | wc -l

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

wget -c http://www.biotrainee.com/jmzeng/tmp/hg38.tss
less hg38.tss
#打开文件后，输入/ID查找
/NM_000546

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

#统计染色体出现的个数则为每个染色体基因的个数
cat  hg38.tss | awk '{print $2}'| sort |uniq -c
#过滤掉碎片基因
cat  hg38.tss | awk '{print $2}'| sort |uniq -c |  grep -v "_"

二十题、解析hg38.tss 文件，统计NM和NR开头的数量，了解NM和NR开头的含义。

cat  hg38.tss | awk '{print $1}' | cut -c 1-2 | sort | uniq -c

ENTRZE ID: 访问NCBI数据库基因的ID号，例如P53基因ID号为：7157
基因NM/NP/NR开头的编号含义 ：基因可以转录为RNA，这个RNA如果可以编码蛋白质则为mRNA,其注释编号为NM_开头，转录的蛋白质注释编号为NP_开头;如果这RNA为非编码RNA（ncRNA）,其注释编号为NR_开头

所以tss文件是啥？