一、在任意文件夹下面创建形如 1/2/3/4/5/6/7/8/9 格式的文件夹系列。
mkdir –p 1/2/3/4/5/6/7/8/9
pwd
二、在创建好的文件夹下面，比如我的是 /tmp/test/1/2/3/4/5/6/7/8/9 ，里面创建文本文件 me.txt。
touch me.txt
三、在文本文件 me.txt 里面输入内容:
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
cat > me.txt
Go to: http://www.biotrainee.com/
I love bioinfomatics.
And you ?
Ctrl+d
四、删除上面创建的文件夹 1/2/3/4/5/6/7/8/9 及文本文件 me.txt
rm –rf 1/
rm –rf me.txt
五、在任意文件夹下面创建 folder1~5这5个文件夹，然后每个文件夹下面继续创建 folder1~5这5个文件夹
mkdir –p folder_{1..5}/folder_{1..5}
六、在第五题创建的每一个文件夹下面都 创建第二题文本文件 me.txt ，内容也要一样。
touch folder_{1..5}/folder_{1..5}/me.txt
vim me.txt (输入内容）
echo folder_{1..5}/ folder_{1..5}|xargs –n 1
echo folder_{1..5}/ folder_{1..5}|xargs –n 1 cp me.txt –v
ls * / * /
七、再次删除掉前面几个步骤建立的文件夹及文件
rm –rf folder_*
rm –rf me.txt
八、下载http://www.biotrainee.com/jmzeng/igv/test.bed 文件，后在里面选择含有 H3K4me3 的那一行是第几行，该文件总共有几行。
wget –c http://www.biotrainee.com/jmzeng/igv/test.bed
ls
pwd
grep -n -o --color H3K4me3 /root/test/test.bed
(-n 标记行数，-o 只显示匹配上的，--color匹配文字出现颜色）
cat test.bed |wc -l
(wc显示文件的行数、单词数、字节数）
九、下载 http://www.biotrainee.com/jmzeng/rmDuplicate.zip 文件，并且解压，查看里面的文件夹结构
wget -c http://www.biotrainee.com/jmzeng/rmDuplicate.zip
ls
yum install -y unzip zip；
unzip rmDuplicate.zip
yum install -y tree
ls
tree rmDuplicate
十、打开第九题解压的文件，进入rmDuplicate/samtools/single 文件夹里面，查看后缀为 .sam 的文件，搞清楚 生物信息学里面的SAM/BAM 定义是什么。
cd rmDuplicate/samtools/single
ls
cat tmp.sam | head
SAM/BAM 定义是什么
Sam:头部注释信息和比对结果，每行包含12列，以tab分隔
Bam：SAM的二进制文件，内存小
1）注释信息可有可无，都是以@开头
@HD，说明符合标准的版本、对比序列的排列顺序；
@SQ，参考序列说明；
@RG，比对上的序列（read）说明；
@PG，使用的程序说明；
@CO，任意的说明信息。
2）详细比对：11个tab隔开的字段
第一列：序列名称
第二列：序列标记信息
1 序列是一对序列中的一个
2 比对结果是一个pair-end比对的末端
4 没有找到位点
8 这个序列是pair中的一个但是没有找到位点
16 在这个比对上的位点，序列与参考序列反向互补
32 这个序列在pair-end中的的mate序列与参考序列反响互补
64 序列是 mate 1
128 序列是 mate 2
如果以上情况都不符合则默认为0
第三列：比对上的参考基因染色体号。没比对上则显示*
第四列：比对上的在参考基因组上的起始位置，没有比对上则为0
第五列：比对质量分数。数字越大，特异性越高。但值为255是指比对质量不可用
第六列：比对结果的CIGAR 的字符串
（即比对的详细情况， 记录插入，删除，错配，后剪切拼接的接头 。以参考序列为基础，使用数字+字母表示比对结果。）
单端测序没有第7-9
第七列：RNEXT ，双端测序中下一个reads比对的参考系列的名称。“＊”是完全没有比对上，“＝”代表完全比对
第八列：PNEXT 如果是双端测序，是指另一端匹配到参考基因组的位置，如果设置为0，那么该列不可用
第九列： TLEN Template的长度，最左边得为正，最右边的为负，中间的不用定义正负，不分区段（single-segment)的比对上，或者不可用时，此处为0；
第十列：SEQ序列片段的序列信息，如果不存储此类信息，此处为’*‘，注意CIGAR中M/I/S/=/X对应数字的和要等于序列长度；
第十一列：QUAL序列的质量信息，格式同FASTQ一样。ASCII码
第十二列：可选字段 格式：TAG:TYPE:VALUE TAG是两个大写字母，每个TAG代表一类信息。TYPE代表TAG对应值的类型（字符串、数组、字节等）
十一、安装 samtools 软件
(提前预装好conda）
conda install samtools
输入`y
十二、打开 后缀为BAM 的文件，找到产生该文件的命令。 提示一下命令是：
/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
cd rmDuplicate/samtools/single
ls
which samtools(查找路径）
samtools(路径) view -H tmp.rmdup.bam
产生bam文件的命令一般都在头部注释信息中，但都会放在最后用
tail
tail -n 3 tmp.header
十三题、根据上面的命令，找到我使用的参考基因组 /home/jianmingzeng/reference/index/bowtie/hg38 具体有多少条染色体。
samtools view -H tmp.rmdup.bam |grep -o -E "SN:chr[0-9]+|SN:chr[a-zA-Z]+"|sort |uniq -c |wc -l
十四题、上面的后缀为BAM 的文件的第二列，只有 0 和 16 两个数字，用 cut/sort/uniq等命令统计它们的个数。
samtools(路径) view tmp.rmdup.bam | cut -f 2 | sort -n | uniq -dc
十五题、重新打开 rmDuplicate/samtools/paired 文件夹下面的后缀为BAM 的文件，再次查看第二列，并且统计
cd rmDuplicate/samtools/paired
ls
samtools(路径) view tmp.rmdup.bam | cut -f 2 | sort -n | uniq -dc
十六题、下载 http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip 文件，并且解压，查看里面的文件夹结构， 这个文件有2.3M，注意留心下载时间及下载速度。
wget -c http://www.biotrainee.com/jmzeng/sickle/sickle-results.zip
unzip sickle-results.zip
tree sickle-results
十七题、解压 sickle-results/single_tmp_fastqc.zip 文件，并且进入解压后的文件夹，找到 fastqc_data.txt 文件，并且搜索该文本文件以 >>开头的有多少行？
unzip sickle-results/single_tmp_fastqc.zip
cd single_tmp_fastqc
cat fastqc_data.txt | grep '>>' | 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
grep NM_000546 hg38.tss
关于它的ID：NM开头的表示标准序列，XM表示预测的蛋白编码序列，NR表示非编码蛋白的mRNA序列，AF开头的表示克隆序列，BC开头的表示模板序列
另外，你可能见过gi|4557284|ref|NM_000646.1|[4557284]这种格式
gi就是代表genebank identifier；ref就是对应的refseq中的ID啦
十九题、解析hg38.tss 文件，统计每条染色体的基因个数。(2种方法)
- less hg38.tss
cut -f2 hg38.tss |cut -d'_' -f1 |sort |uniq -c |sort -rn - cat hg38.tss |cut -f 2|grep -o -E "chr[0-9]{1,2}|chr[a-zA-Z]{1,2}"|sort |uniq -dc
二十题、解析hg38.tss 文件，统计NM和NR开头的序列，了解NM和NR开头的含义。(3种方法）
- cat hg38.tss |grep 'NM' | wc -l
cat hg38.tss |grep 'NR' | wc -l - grep -o -E "NM|NR" hg38.tss |sort |uniq -dc
-
cat hg38.tss |awk '{print $1}'|cut -c1-2|sort|uniq -c
(只带有前2个字母进行统计）
NM:开头表示标准序列。可以转录成蛋白质的基因
NR:开头表示非编码蛋白的mRNA序列
数据练习
- (首先下载好数据）
- 看文件大小
ll -h *.gff - 看文本内容
less -SN *.gff - 去除多余的#注释行与 空行，并打印行号
grep -v "#" *.gff | grep -v "^$" | wc -l - 截取文件1-5列，将第2列去除，输出默认的前10行到test.txt中
cut -f 1,3,4,5 *.gff > test.txt
cat test.txt | head
- (然后对截取的test.txt进行处理）
- 根据第2列的feature进行排序
sort -k 2,2 test.txt | head -n5 - 先根据第1列Chr数字大小降序排序(第一个字段的第四个字符），再根据第2列排序
sort -k 1.4,1.4nr -k 2,2 *.txt | head -n5 - 统计整体的feature(第3列）
cut -f 3 *.gff |sort|uniq -c >feature.txt - 哪个feature最多
sort -k 1r feature.txt
练习：
GENCODE下载人类基因组GRCh38注释gff3，然后统计人类基因组feature（第1列 参照序列，2来源，3类型，4起点，5终点，6得分，7链（正负链），8步进，9属性）
1.下载
wget -c ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_28/gencode.v28.annotation.gff3.gz
2.解压
gunzip gencode.v28.annotation.gff3.gz
3.看文本内容
less -SN *.gff3
4.去除多余的#注释行与 空行，统计整体的feature(第3列）
grep -v "#" *.gff3 | grep -v "^$" | cut -f 3 |sort|uniq -c >feature.txt
5.查看文件(前10行）
cat feature.txt | head
- 哪个feature最少（从小到大）
sort -k 1n feature.txt
(从大到小）
sort -k1,1nr feature.txt
